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INTRODOCTION 


E1CKGECDNE 


The  Fcurier  transform  has  been  with  us  many  years  and 
the  uniqueness  of  the  transform  allows  Fourier  analysis 
technigues  developed  in  one  field  to  be  applied  to  many 
diverse  areas.  lypical  application  areas  of  the  Fcurier 
transform  include: 

Linear  Systems  -  The  Fourier  transform  of  the 
output  cf  a  linear  system  is  given  by  the  product 
cf  the  system  transfer  function  and  the  Fourier 
transform  of  the  input  signal. 

Antennas  —  The  field  pattern  cf  an  antenna  is 
given  by  the  Fourier  transform  of  the  antenna 
current  illumination. 

Optics  -  Optical  systems  have  the  property  that  a 
Fourier  transform  relation  exists  between  the 
light  amplitude  distribution  at  the  front  and  back 
focal  planes  cf  a  converging  lens. 

Randcm  Process  —  The  power  density  spectrum  of  a 
random  process  is  given  by  the  Fourier  transform 
cf  the  auto-ccrrelation  function  cf  the  process. 

Probability  —  the  characteristic  furction  of  a 
randcm  variable  is  defined  as  the  Fourier 
transform  of  the  probability  density  function  of 
the  randcm  variable. 

Quantum  Physics  —  The  "uncertainty  principle"  in 
quantum  theory  is  fundamentally  associated  with 
the  Fourier  transform  since  particle  momentum  and 
position  are  essentially  related  through  the 
Fourier  transform. 

Eouncar j-Value  Problems  -  The  solution  of  partial 
differential  eguations  can  re  obtained  by  means  of 
the  Fctrier  transform. 
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Although  Fourier  analysis  allows  one  to  more  easily 
examine  a  functicn  from  another  point  of  view, i.e.,  the 
transform  dciain,  the  methods  available  before  1S65  to 
compute  the  Pcurier  coefficients  were  very  time  consuming 
and  costly.  Then  in  1965  Cocley  and  Tukey  published  their 
mathematical  algorithm  which  computed  the  Fcurier 
coexficents  with  much  less  computation  effort  than  had  teen 
required  in  the  past.  This  method  has  become  known  as  the 
"fast  Fcurier  transform"  or  FFT  and  has  produced  major 
changes  in  computational  techniques  used  in  digital  spectral 
analysis,  filter  simulation  and  related  fields. 

Although  the  FFT  is  the  most  well  known  algorithm  to 
compute  the  fcurier  coefficients,  it  is  net  the  only  method 
fcr  computing  the  "discrete  Fcurier  transform"  (DFT)  .  The 
chirp-z-transf crm  (CZT)  is  an  algorithm  fcr  evaluating  the 
z-transfcrm  cf  a  finite  duration  sequence  along  certain 
general  ccntcurs  in  the  z-plane.  Evaluating  the  DFT  cf  the 
seguence  is  a  special  case  of  the  CZT,  i.e.,  the  z-transform 
cf  the  finite  duration  seguence  is  evaluated  on  the  unit 
circle  in  the  z-plane.  Although  the  CZT  is  not  guite  as 
efficient  as  the  FJT,  it  eliminates  many  of  the  restrictions 
of  the  latter.  These  restrictions  are  listed  at  the  end  of 
the  CZT  algorithm  derivation.  The  "prime  transform"  is  an 
algorithm  that  fellows  directly  from  the  CZT  derivation.  It 
is  based  upon  "N",  the  number  of  samples  cf  data,  being  an 
edd  prime. 

The  CZT  and  prime  transform  algorithms  are  receiving 
consideratle  attention  for  applications  in  spectral 
analysis.  lhis  has  been  brought  about  because  the  bulk  of 
the  computations  of  these  algorithms  is  performed  ty  a 
transversal  filter,  which  is  easily  implemented  hy  charge 
transfer  devices  (CTD)  or  surface  acoustic  wave  (SAW) 
devices.   Hardware   and  software  implementations  of  the  FFT, 
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CZT,  and  prime  transform  algorithms  are  discussed   later   in 
this  chapter. 


EASIC  CCKCEPIS  OF  FCORIER  TRANSFORM  ANALYSIS 


The  esseice  of  the  Fourier  transform  of  a  waveform  is  to 
decompose  or  separate  the  waveform  into  a  sum  of  sinusoids 
cf  different  frequencies.  The  pictorial  representation  of 
the  Fourier  transform  is  a  diagram  which  displavs  the 
amplitude  and  frequency  of  each  of  the  determined  sinusoids. 

Mathematically,  the  AMPLITUDE  SPECTRAL  DENSITY  or  more 
generally  the  fourier  transform  of  x(t)  is  given  b$  the 
relationship 


-j2TTft 

X  (f)  =   x(t)  e       dt  (1-1) 


where  x  (t)  is  the  waveform  to  be  decomposed  into  a  sum  of 
sinusoids,  X  (f )  is  the  Fourier  transform  of  x  (t)  ,  and  j=y-T. 
The  INVEESE  fCURIEE  TRANSFORM  is  given  by  the  relationship 


/        +j2TTft  „  „ 

x(t)  =    X.(f)e      df  (1-2> 


The   Fouiier   transform  and  its  inverse  may  conveniently  be 
written  as 

v(t)  <3zz£>  v(f) 
FOURIER  TEANSFCRM  EAIR 
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For  discrete  signal  processing  the  equations  take  en  the 
following  form.  For  a  periodic  function  x(t)  the  Fourier 
tiansfcra  is  defined  by 


N-1 

"\      -j(2TT/N)nk 
X(k)  =  >x(n)e  ,  k=0,1,...,N-1    (1-3) 

n  =  0 


The  "true"  freguency  component  is  given  by 


N-1 

V~     -j(2TT/N)nk 
X(k/NT)  =  Nx(nT)e  (1-4) 

n=0 


Egs.[1-3]  and  [1-4]  are  known  as  the  EISCRETE  FCUEIEB 
TEiKSFOBH  (£11).  The  inverse  DFT  is  given  by  the  following 
relation 


N-1 

V~  +j(2TT/N)nk 

x(n)    =    (1/N)  ^X(k)e  ,n=0,1,...,N-1  (1-5) 

k  =  0 


and    for    "true"   time 


1-1 


♦J(2tt/N)  nk 
x(nT)    =     (1/N))     X(k/NT)e  (1-6) 


k=0 
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where 


N  — >  total  number  of  samples  cf  time 
and  samples  of  frequency 

T  — >  sample  period 

n  — >   sample  period  index 

k  — ^  harmonic  index 

P  =  (1/NT)  ->  fundamental  frequency 


Fcr  future  reference,  the  basic  properties  cf  the 
continuous  Fcurier  transform  and  discrete  Fourier  tracsform 
are  sumiarized  in  Iable-I.  Fig.[1]  gives  a  pictorial 
representation  of  some  of  the  most  commonly  used  Fcurier 
transfer!  pairs. 
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IKFLFKFMMIGN  PHCELESS 


lh€   two   problems   most   often  encountered  in  using  the 
discrete  Jcirier  transform  are  aliasing  and  leakage. 


1 .   Jliasi.ng 


Aliasing  refers  to  the  fact  that  the  high-frequency 
components  cf  a  time  function  can  impersonate  low 
frequencies  if  the  sampling  rate  is  too  low.  This  problem 
is  corrected  by  sampling  the  signal  at  a  rate  >  the  Njguist 
Freguency,  i.e.,  demanding  that  the  sampling  rate  be  high 
enough  for  the  highest  frequency  present  to  be  sampled  at 
least  twice  during  the  period.  Fig.[2]  gives  a  pictorial 
representation  of  a  signal  without  aliasing. 
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figure  2  -   Fourier  Transform  of  a  Wavefcrn  Sampled  at  the 
Nyguist  Sampling  Rate. 


Iten  the  sampling  period  T  of  Fig.[2-c]  is  increased 
as  shown  id  fig.[3-c:],  the  equidistant  impulses  of  (f) 
become  acre  closely  spaced,  Fig. [3-d].  Eecause  of  the 
decreased  spacing  of  the  freguency  iapulses,  their 
convolution  with  the  freguency  function,  a (f ) ,  Fig.[3-b 
results  in  tke  overlapping  waveform  illustrated  results  in 
the  overlapping  waveform  illustrated  in  Fig.[3-f].  This 
overlap  is  the  result  of  aliasing  by  the  high  freguency 
ccupcnents. 

Ike  aliasing  phenomenon  will  occur  anytime  a  sampled 
analog   signal   contains   freguency  components  greater  than 

half  the  santling  freguency,  (f  ) /2.   The  aliasing  frecuency 

s 

relationship   is   such   that  all  frequencies  higher  than  the 

analysis  bandwidth  [f    =  (f  )/2]   "fold"   dcwn   into   that 

max     s 

band   in   an   accordian-like   manner   with   a   fold  at  every 
multiple  cf  f 
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Than  the  Nyguist  Sampling  Bate. 
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2.   Ieakace 

If  a  periodic  ,  band  limited  function  is  sampled  and 
truncated  tc  consist  of  other  than  an  integer  multiple  of 
the  period,  the  resulting  discrete  and  continuous  Fourier 
transforms  %ill  differ  considerably.  This  problem  is 
inherent  in  the  Fourier  analysis  of  any  finite  record  of 
data.  The  record  has  been  formed  by  looking  at  the  actual 
signal  fcr  1  seconds  and  by  neglecting  everything  that 
happened  hefcre  cr  after  this  period.  As  shewn  in 
Fig.[4-cj,  this  is  eguivalent  to  multiplying  the  signal  by  a 
rectangular  window. 

If  the  continuous  Fourier  transform  of  the  pure 
cosine  wave  cf  Fig.[4-a]  had  been  found,  its  contritution 
would  ha^e  been  limited  to  a  pair  of  impulse  functions  en 
the  frequency  axis,  Fig.[4-a]. 
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Ihe  effect  cf  truncation  at  other  than  a  multiple  of 
the  pericd  cf  tte  sampled  signal  is  to  create  a  periodic 
function  *ith  sharp  discontinuities  as  shewn  in  Fig. [5], 
The  multiplication  by  the  data  window  in  the  time  domain  is 
eguivalent  tc  performing  convolution  in  the  freguency 
domain.  Ccnseguen tly ,  the  freguency  function  is  no  longer  a 
single  inpulse  but  rather  a  continuous  function  of  freguency 
with  a  lecal  maxiium  centered  at  the  original  impulse  and  a 
series  cf  spurious  peaks  termed  sidelobes.  These  sidelcbes 
are  responsible  for  the  additional  freguency  components 
"leakage*',  which  cccur  after  freguency  domain  sampling.  An 
expanded  view  cf  sampling  in  the  freguency  domain  is  shown 
in  fig. [6], 
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figure  5  -   EFT  cf  a  Periodic  Waveform: Truncation  Interval 
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Figure  6  -   Expanded  Illustration  cf  the  Convolution  of 
Figure  [ 5e]. 


Ihe  tsual  approach  to  reduce  the  ancunt  of  "leakage" 
through  ttese  sidelotes  consists  cf  applying  a  data  window 
or  weighting  function  to  the  time  series,  which  has  lower 
sidelobes  in  the  frequency  domain  than  the  rectangular  data 
window.  The  process  of  weighting  consists  of 
affitlitude-ncculatirg  the  signal  prior  to  Fourier  analysis. 
The  weighting  process  is  illustrated  in  Fig. [7], 
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FOURIER  TRANSFORM 
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TRANSFORM  OF  THE 
TIME  FUNCTION 
CONVOLVED  WITH 

TRANSFORM  OF  THE 
WEIGHTING  FUNCTION 


Ficure  7  -   The  Weighting  Function  Frocess 
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Ihere  are  many  such  windows  but  one  of  the  sinplist 
and  most  often  used  is  the  Hanning  function.  This  function 
is  a  ccsine-bell  curve  given  by 


w  (t)  =  1  -  1ccs2TTt 
I   2   ~Tc 


0<t<Tc 


(1-7) 


where  I   is  the  truncation  interval.   The  magnitude   of   the 
c 

fcurier  transform  cf  the  Hanning  function  is  given  by 


|i(i)  I  =  JQ(f)  *  JCQU+J  )  +  C(f-1  )  ] 
^      a     Tc       Tc 


(1-8) 


where  Q(f)  =  [  sin  Cttc£  )  ]/TTf 

The   Hannirg   function   Fourier   transform   fair  is  shewn  in 
Fie. [8]. 

Windowing  of  the  input  frame  by  the  Banning  function 
reduces  "lea>age"  ty  imposing  a  guasi-pericdicity  en  the 
input  sicral  and  eliminating  the  discontinuities  between 
periods.  The  Hanning  function  widens  the  nain  lobe  of  the 
(sin  x)/x  shape  and  reduces  the  spectral  amplitude  ty  a 
little  mere  than  4dE,  but  it  causes  the  "interfering" 
sidelcbes  to  fall  off  at  18d3  per  octave  instead  cf  6dE  as 
in  the  rectargular  window. 
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It  shculd  be  remembered  that  the  non-zero  frequency 
ccuponents  are  considerably  EROADENED  cr  SMEABEE  with 
respect  tc  tfce  iapulse  function.  In  general,  the  mere  the 
leakage  is  reduced,  the  broader  or  more  smeared  the  results 
of  the  discrete  Fourier  transform  appear.  Iable-II  contains 
scae  of  the  Here  cemmenly  used  data  windows. 
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TABLE-II 
Data  Windows 


WEIGHTING  FUNCTION 


RECTANGULAR 


1/B 


TRIANGLE 


3-dB 

BW 


0.85B 


1.253 


NOISE 
BW 


IB 


1.35B 


HIGHEST 
SIDELOBE 
LEVEL  (dB) 


•13 


■26 


ASYMPTOTIC 
ROLL-OFF 
(dB /OCTAVE) 


12 


COSINE 


1.17B 


1.26B 


•23 


12 


HANNING  =   (COSINE)2 


l.UB 


1.5B 


■32 


18 


(cosine) 3 


1.6IB 


1.73B 


■39 


24 


1.88B 


1.9B 


-48 
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HAMMING  =  (COSINE)2  + 

3#  PEDESTAL 


1.3B 


1.36B 


-42 


6d£/CCT 

BEYOND 

5B 
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C.   SFECTBUE  iNAIYZEB  PERFORMANCE  FIGURES  OF  MERIT 


1  •   Saj^l inq  Frequency  and  Antialiasing  Filter 

The  sampling  frequency  must  satisfy  the  Nyguist 
sampling  ttecrem  or  aliasing  will  result.  Ecr  data  that  has 
conpcnents  bejcnd  the  analysis  range,  the  components  are 
attenuated  ty  means  of  a  lowpass  filter,  known  as  an 
"anti-aliasing"  filter.  A  typical  "fall  off"  rate  for 
anti-aliasing  filters  in  today's  equipment  is  120  dE/octave. 
In  practice,  a  sampling  rate  of  approximately  three  times 
the  width  cf  the  desired  analysis  range  is  a  good 
ccmpromise.  Sampling  at  higher  rates  requires  more 
equipment,  *hich  is  undersiratle  for  eccncmic  reascns. 
Sampling  at  lower  rates  fails,  to  take  into  account  the 
actual  fall-cff  of  available  anti-aliasing  filters,  anc  thus 
dees  not  eliminate  aliasing  completely. 

2«   l^ad width  (Frequency  Ban^e). 

Tie  bandwidth  determining  parameters  of  a  spectrum 
analyzer  are: 

1.  lie  memory  capacity  (which  determines  the 
numter  of  samples  available  for  processing) . 

2.  The  samcling  frequency  (which,  in  combination 
*itk  1',  determines  the  length  cf  signal  (ie 
seconds)  which  is  available  fcr  processing) . 

~.  The  shape  of  the  "windew  function"  also  known 
as  the  "weighting  function". 
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liiCluticn  cr  3-dB  Bandwidth 


In  Table-II  the  3-dB  bandwidth  for  various  weighting 
furcticns  were  compared  using  the  relation 


B(H2)  = 


Seighting-functicn  length  (sec) 


The  3-dE  tandwidth  is  a  good  measure  cf  a  spectrum 
analyzer's  ability  to  resolve  two  equal-amplitude  adjacent 
sire  waves.  Studying  Table-II,  it  appears  that  the 
rectangular  weighting  function  has  the  best  resoluticr,  but 
in  reality,  the  response  characteristic  is  limited  in 
usefulness  ty  the  side-lobe  structure.  In  practice,  it  is 
mere  important  to  resolve  unegual  than  equal  amplitude 
adjacent  frequency  components. 

4»   Seise  Bandwidth 

Scne-times  called  "effective  bandwidth",  it  is  used 
to  convert  cr  normalize  power  spectrum  measurements  to  power 
spectral  density  (power  per  unit  frequency) .  The  noise 
bandwidth  is  the  bandwidtn  cf  a  hypothetical  rectangular 
filter,  fchich  passes  a  signal  with  the  same  mean-square 
value  as  the  actual  filter  when  the  filter  input  is  white 
Gaussian  rcise. 

5»   Eia^iSi  Sidslobe  Level  an.d  Asymptotic  Rcll-of f 

It  is  useful  to  specify  these  parameter,  since  it 
provides  seme  indication  of  the  prominence  cf  the  side-lobe 
structure.   The  largest  sidelobe  is  generally  the  first  cne, 
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with  the  Banning  function  being  amcng  the   exceptions.    The 
third  side  lcte  of  the  Hamming  function  is  the  largest. 

Ihe  presence  of  side  lores  makes  it  difficult  to 
distinguish  between  a  high-level  frequency  component  and  a 
lew-level  component  which  is  close  in  frequency.  For  this 
reason,  rectangular  weighting,  which  leads  to  a  spectrum 
response  egual  to  the  true  Fourier  transform  of  the 
tine-truncated  signal,  is  used  only  occasionally  when 
analyzing  transients  or  shocks. 


6 .   2i5SlizSided  Spectra 

Although  it  is  net  an  actual  measure  of  performace, 
a  discussion  of  single-sided  spectra  is  warranted  at  this 
tine  so  that  the  results  obtained  frcm  the  various 
algorithms  is  fully  understood. 

The  finite  discrete  transform  presumes  that  the 
input  waveform  is  the  sum  cf  cosine-phase  and  sine-phase 
freguency  ecapenents,  all  harmonically-related  to  the  lowest 
frequency  that  can  be  recognized.  That  fundamental  is  the 
freguency  that  completes  one  cycle  in  the  length  of  the 
input  frame.  The  calculated  harmonics  range  frcm  <d-c) 
through  1  (the  fundamental) , 2,3, etc. ,  to  N/2,  the  maximum 
recognizable  harmonic,  which  completes  one  cycle  in  the  two 
intervals  associated  with  two  pents  of  the  input  (i.e.,  two 
San  pies  per  cycle)  . 

In  the  general  case,  the  discrete  Fourier  transform 
can  operate  en  a  complex  time  series  input  function  and 
produce  a  cenplex  spectrum.  In  such  a  case,  an  input 
function  cf  B  complex  values  (i.e.,  2N  independent  real  and 
imaginary  parts)  yields  a  spectrum  of  2N  independent  real 
and    iiaginary    parts  representing   N   complex   frequency 
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coefficients  centered  about  0  (d-c)  and  ranging  from  harmonic 
-N/2  to  +  N/2.  The  "negative-frequency"  components  are  an 
integral  and  necessary  part  of  the  function. 

A  real  input  function,  as  an  electrical  waveform, 
may  te  considered  for  this  general  case  as  a  complex 
futcticn  whose  imaginary  parts  are  all  zero.  The  discrete 
transfers  of  such  an  input  yields  a  complex  spectrum  as 
afceve.  For  this  special-case  input  however,  the  spectrum 
will  be  symmetrical  about  the  center  or  d-c  point.  The 
necative-f recuency  parts  from  0  to  -N/2  are  an  image  of  the 
positive  f reguencies-,  the  real  or  cosine  parts  having  even 
symmetry  arc  the  imaginary,  sine  parts  odd.  The 
necative-f recuency  half  thus  becomes  redundant,  and  an  input 
of  N  real  points  yields  N  independent  real  and  imaginary 
parts  of  N/2  complex  spectral  coefficients. 

Nonetheless,  the  negative-frequency  coefficients  are 

still   a   necessary   part   of   the   spectrum   in   that   they 

represent  half  the  amplitude  of  the  harmonic  components.   In 

ether  words,  if  the  input  function  contains   a   cosine-phase 

component    cf   amplitude   A   at   frequency   3f  ,   ther   the 

c  0 

resultant  positive-frequency  coefficient  for   frequency   3f 

will    have    amplitude   A/2    and   the   negative-f recuency 

th 
coefficient  fcr  -3f   will  also  have  amplitude  A/2.   The   0 
0 

or   d-c   coefficient,   however,   dividing   or  being  "shared" 

between  the  pesitive  and  negative  halves  cf  the  spectrum, 
will  have  full  amplitude  A  for  an  input  d-c  level  of  A, 
which  is  the  pcint  cf  this  discussion. 


36 


SPEC1SAL  ANALYSIS  ALGORITHMS 


1  •  lk£   I^st  £eurier  Transform 

lte  fast  Fourier  transform  (FFT)  is  a  highly 
efficient  procedure  fcr  computing  the  DFT  cf  a  time  scries. 
It  takes  advantage  cf  the  fact  that  the  calculation  of  the 
ccefficierts  cf  the  DFT  can  be  carried  cut  iteratively, 
which  results  in  a  considerable  savings  of  computation  time. 
Tc  gain  icsicht  to  the  FFT  algorithm,  a  general  developement 
is  preserted  followed  by  several  examples.  A  more  in  depth 
study   of  the  Ell  may  be  found  in  Eefs.[6  ],  [16],  and  [17]. 

a.   General  Cevelopement  cf  the  FFT  Algorithm 

lie  discrete  Fourier  transform  Eg. [1-3]  is 


,  k=0,1,...,N-1    (1-3) 


which  ma}'  be  written  as 


(1-9) 
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wtere 


H  =  € 


J2TT/N 


(1-10) 


nk 
is  can  be  seen,  W    is  periodic  and  of  period  N: 


1.6 


(D  +  mN)  (k+lN)     nk 

H  =  1 


,1  =  0,±1,... 


(1-11) 


Id   carrying   cut   the  FFT   algorithm,   the   symmetry    and 

nk 
periodicity   cf  n        are  exploited  to  achieve  an  increase  in 

efficiency.   Tc   emphasize   this   periodicity,   W   will   be 

rewritten  as  i   in  carrying  cat  the  general  developemert . 

nk 
In  Eg.  £1-9],  if        is  a  complex  exponential,   thus 

an  examination  of  Ig-CI-9]  shews  that,  in  the  case  when  x  (n) 
is  a  ccmclex  sequence,  a  complete  direct  evaluation  cf  an 
N-pcint  EFT  requires  (N-1)2  complex  multiplications  and 
N<H-1)  cciflex  additions.1  Fig. [9]  shews  the  compariscn  of 
multiplications  required  by  the  direct  calculation  cf  the 
Ell  and  fcr  the  base  2  FFT  algorithm. 


*In  most  literature,  N  is  assumed  to  be   large, 
sc  that  (N-1)  can  te  approximated  by  N. 
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Figure    9    -      Multiplications    Required   to    Compute    the    BIT, 
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t.   Eecimation-In-Time  (DIT)  Algorithm 

Ibe  general  principle  behind  the  PFT  is  to  break 
the  original  N-pcint  sequence  into  twc  shcrter  seguences, 
the  DFT's  c£  which  can  be  combined  to  give  the  DPT  of  the 
original  N-pcint  seguence. 

Assuming   N   is   a   power   of   2,    define    two 

(N/2) -point   segcences   x  (n)   and  x  (n)  as  the  even  and  odd 

1  2 

meabers  cf  the  seguence  x(n)  respectively,  i.e., 


x1(n)  =  x(2n) 
x«(n)  =  x(2n+l 


N 


N 


-  1 


(1-12) 


n  =  0,1, ...,§  -  1 


The  N-point  DFT  of  the  sequence  x(n)  can  be  written  as 
N-l  N-l 

X( 


[k)  =  Yx(n)wjf  +  Vx(n)wf 


n=0 
n  even 


n=0 
n  odd 


(1-13) 


N/2-1 


N/2-1 


>x(2n)W 


.2nk 
N 


^x(2n+l)wi2n+1)k 


(1-1*0 


Then  using  the  following  identity 


W2 
WN 


J(2tt/N) 


2  _  J 


=  e 


[2<n/(N/2 


=  W 


N/2 


Eq. [l-l^l  can  written  in  the  form 
N/2-1  N/2-1 

:1 


IWN/2   +  /   x2(n)W 


nk 
N/2 


n=0 


(1-15) 


(1-16 
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X(k)  =  x1(k)  +  w*x2(k) 


(1-17) 


Each   cf   the  sums  in  Eg. [1-16]   is  an  N/2-pcint 

Efl,   the   first   sujd   being   the   N/2   point   DFT   cf    the 

even-numterec   points   of   x  (n)  and  the  second  being  the  N/2 

ccint   Ell   cf   the   odd-numbered   points   cf   the   original 

seguence.    ilthcugh   the   index   k   ranges   over   N  values, 

k=C,  1 , . .  .  ,N-  1 ,  each  of  the  sums  need  only  be  computed  for   k 

between  C  ace  N/2-1,  since  X  (k)  and  X  (k)  are  each  periodic 

1         2 

in  k  with  period  N/2.  After  the  two  DPT*s   corresponding   to 

the   two   suns   in   Eg. [1-16]   are   computed,   they  are  then 

ccabined   to   jield   the   N-point   DFT,   X  (k)   as   given   by 

Eg. [1-17]. 
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2-   lie  Cfcj.r£-Z-Transforffl  JCZJ1 

The  chirp-z-transform  (CZT)  is  an  algorithm  for 
evaluating  th€  z-transform  of  a  finite  duration  sequence 
alcng  certain  general  contours  in  the  z-plane.  The 
following  derivation  is  based  upon  a  reccrt  by  Bariner, 
Schafer,aEd  Eader  [18]. 

a.   lie  Derivation  of  the  CZT  Algorithm 


In   gereral,   the   z-transfcrm   cf  a  sequence  of 
sanples  |x(n)}  is  defined  as 


(n)z  Q  (1-18) 


Assumminc   ttat   the   right  side  of  Eg. [1-18]   converges  for 

some  values  cf  z,  and   restricting  the   evaluation   of  the 

z-transfcrn    to    a    finite  duration   sequence  of 
sanples,     Ig.[1-18]. 


N-1 
X  <2)  =  )   x(n)z  (1-19) 

n  =  0 


The  special  case  cf  the  z-transform  which  has 
received  considerable  attention  is  the  set  cf  points  equally 
spaced  arcund  the  unit  circle  of  the  z-plane,  i.e.r 
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Z   =  exp  (j2TTk/N)  ,    k=0,1,...,N-1 

K 


(1-20) 


which  gives 


N-1 
X(2)  =  \  x(r 

n  =  0 


-  j  (2tt/N)  nk 


,  k=0,1,...,N-1   (1-21) 


By  comparing  Eg. [1-2.1]  and  Eq.[1-3],  cne  can  see  that  it  is 
th€  discrete  Fcuiier  transform  (DPI) .  The  DPI  can  be 
evaluated  mere  efficiently  by  using  the  FFT  rather  thar  the 
CZT,  tut  the  latter  allows  cne  to  evaluate  Eq.x-19]  along 
the  mere  cereral  ccntcur 


-k 


Z   =  AW 
k 


,  k=o,i,...,a-i 


(1-22) 


where  a  is  at  arbitrary  integer  (not  necessarily  equal  to  N) 
and  A  and  *  are  arbitrary  complex  numbers  defined  by 


A  =  A    exp  ( J2TT9  ) 
c        o 


1-23) 


fc  =  »  exp  (J2tt^  ) 
o        o 


(1-24) 


Fig. [101   shews  the  significance  of  A  ,   W  ,  9  ,   and  0       in 

a  L   J  o    c    c         0 

the  z-plare. 
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Z -PLANE 


figure  10  -   Contour  of  the  CZT  in  the  Z-Plane 


U4 


fcr  the  special  case  of  the  Dfl,  A=1  ,  M=S,  and 
H=exp  (-j2tt/N)  ;  thus  the  DPT  coefficients  of  a  finite 
duration  secuence  are  the  values  of  the  z-transform  of  that 
sane  sequence  at  N  evenly  spaced  points  around  the  unit 
circle. 


lhe  parameter  W   determines  the   rate   at   which 
o 

the   contctr   spirals:   if   W   is  greater  than  unity,  the 

o 

contour  spirals  toward  the  origin  as  K  increases,  and  if   tf 

o 

is   less   than   unity,   the   contour   spirals   outward   as  K 

increases.   lhe  parameters  A   and  2tt9  are  the  locaticn   in 

o        o 

radius   and   angle,  respectively,  cf  the  first  sample,  i.e., 


fcr  k=0.   lie  remaining  samples  are  located  alcng  the  spiral 

ccntcur   *itt   an   angular   spacing  of  2ttJ2'  .   The  equivalent 

c 

s-tlane  ccntcur  of  Fig. [10]  is  shown  in  Fig. [11]. 

S- PLANE 


AO  |  in  WQ  -H  ^ 


Co  ~  T      0 


Figure  11  - 


o    T" 
Contour  cf  the  CZT  in  the  S-Plane 
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lfc€  equivalent  s-plane  contour   beginning   point 
may  be  found  by  iraking  the  substitution 


th€D 


cr 


sT 

Z  =  € 


z   =  A  =  exp  (s  T) 
c  o 


s   =  (1/1)  InA 
c 


=    a      ♦  jw   =  (1/T)  (InA   +  J2u9)  (1-25) 

CO  0  0 


For  a  genera]  point  en  the  s-plane  contour 

-k 


tt€D 


Z   =  exp  (S  T)  =  AH 

k        x 


-k 
s   =  (1/T)  la  (AW   ) 
k 


anc  finally 

s   =  (1/1)  (InA-klnW)  =  s  -(k/T)ln« 
k  c 

k=0,1,...,M-1       (1-26) 

Ihcs   the   spiialing   contours   in  the  z-plane  correspond  to 
straight  lines  in  the  s-plane.   The  rate  cf  spiraling  in  the 
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z-tlane  deteraines  the  slope  cf  the  lire  in  the  s-plane. 

Ihe  task  now  at  hand   is   to   compute   Eg. [1-19] 

along   a  general  contour  given  by  Eg. [1-22],   The  evaluation 

of  X(z)  at  z  =  z  ,  will  be   written  using   the   notation   X  , 
k  k 


1-1 


x  (n)  z 


H-27) 


n=0 


or 


N-1 

\  '    ~n  nk 
X   =  \x(n>A   W      ,  k=0,1,  ...  ,K-1  (1-28) 


n=0 


where  N  is  tie  length  of  the  sequence  x(n). 

taking  the   following   substitution,   which   was 
crigninated  rj  Eluestein  [ 5 ] - 


nk  =  na  ♦  k2  -  tk-n)  2 
2 


(1-29) 


gives 
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N-1 

y-n   [n2  +  k2-  (k-n)  ]/2 
x<n)A      WL  J/ 


n=0 


1-1 


■2 

n=0 


-n  n2/2   k2/2  -  (k-n) 2/2 
[s(n)A   H     ]W     I  (1-30) 


although  Eg. [1-30]  appears  to  be  mere  complicated  than 
Eg. [1-28],  it  allows  efficient  hardware  inplementation. 
Eg. [1-30]  nay  te  viewed  as  basically  a  three  step  process: 


1.   lie   first    step    consists    of    a    complex 
pre-iultip lication  forming  a  new  sequence 


j  (C)  =  x(n>A  Dwn   "  ,n=0,1,. .. ,N-1  (1-31) 


2 


Ibe   seguence   y(n)   is   then   convolved   with   the 
seqc€nce  v(n)  derined  as 


-n2/2 
v  (E)  =  N  (1-32) 


tc  give  the  seguence  g (k)  defined  as 

N-1 
g(k)  =    y(n)v(k-n)    ,k=0, 1  , . .  .  ,fl-1         (1-33) 

n  =  0 

k2/2 

3.   The  seguence  g(k)  is  then  post-multiplied  fcy   W 
tc  cive 
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X  =  w 
k 


k2/2 


N-1 

\y(n)v(k-n)  ,k=0,  1 , . . .  ,H-1 

n=0 


H-34) 


A   pictorial  representation  of  the  prccessiDg  operations  of 
tb€  CZI  algorithm  is  shewn  in  Pig. [12]. 


A~nWn  /2 


figure  12  -   Processing  Operations  of  the  CZT  Algorithm, 
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The  advantages  cf  the  CZT  algorithm  ever   the   standard 
PIT  are  suanarized  below: 

1.  K,  the  number  of  points  in  the  input  sequence, 
dcesn*t  have  to  be  egual  to  M,  the  number  of 
points  ct  which  the  z-transform  is  evaluated. 

2.  Neither  N  nor  M  has  tc  be  composite  -  both  can 
be  prime  numbers. 


3.  The  starting  point  in  the  z-plane  and  the 
angular  spacicg  or  tne  z  »s  are   arbitrary.    Thus 

k 
the    freguency   resolution   and  the  range   of 

fregcencies  are  arbitrary. 

4.  lie  contour  does  not  have  to  be  a  circle  in 
the  2-plane.  The  CZT  algorithm  can  evaluate  the 
z-transfcrm  along  spiral  contours  inside  and 
outside  the  unit  circle  of  the  z-plane.  By 
evaluating  the  z-transform  along  a  selected 
contour,  the  transfer  function  of  a  linear  systen 
can  te  adjusted  to  sharpen  or  flatten  the 
freguency  response.  The  choice  of  contour  is 
dependert  cd  which  of  the  system's  poles  should  be 
emphasized. 

5.  If  1*1,  I*-Hf  W=  exp  (-i2iT/N)  ,  then  the  CZT  can 
be  used  tc  evaluate  the  DFT,  even  when  N  is  prime. 

6.  Kith  the  CZT  algorithm^  interpolation  between 
time  saifles  cf  a  Bandliaited  functicn  can  be 
obtained  with  greater  ease  with  the  DFT,  and 
interpolation  at  arbitrary  sampling  intervals  can 
be  performed  efficiently. 
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3.   The  Erisae  Transform 


a.   Tte  Prime  Transform  Derivation 

The  fcllpwing  derivation  is  based  on  a  paper  by 
Shitehouse,  Ceans,  and  speiser  [26],  A  brief  review  cf  the 
nuater  tbeorj  properties  used  in  this  derivation  is  given  in 
Appendix  B. 

In  the  specific  case  when  H  is  an  odd  prime. 
Eg. [1-3]  can  te  written  as 


N-1 

Vx(n) 


(1-35) 


2 


n  =  0 


tc    give    the    (D-C)     freguency    component   and 


N-1 

Y-j2TTnk/N 
x(n)e  ,k=1,...,N-1  (1-36) 


n=1 


tc  give  tte  remaining  components. 

Since  cnly  non-zero  values  of  n  and  k  occur  in 
tte  right  bacd  side  of  Eq.£1-36],  it  is  possible  tc  replace 
the  product  tk   by  a  primitive  root  raised  tc  a  sum  [19]. 

Proceeding  ,  in  the  summation  cf  Eg. [1-36],  the 
terms  are  permuted  and  the  order  of  the  equation  changed   by 
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the    fcllciiiDc    transformations: 

(n)p  —  ((r(n>P))  =  rn  modulo  N 
(k)p  —  ((r(k)P))  =  rk  modulo  N 


(1-37) 


n,k=l,...,N-l 
where  "r"  is  a   primitive   root   and   "p"   implies  permuted. 

N-1  0 

Ncting      that      r  modulo      N      =      r        modulo  Nf    Eg. [1-36]   is 


viritten    as 


N-1 
\x(n)         ,„,„      exp(-^r((n^(k)P) 


((r(k)P)) 


((r(n)P)) 

n=0 


(1-38) 


An  examination  of  Eq.i  1-38 1,   shows  that  the  sequence 

/v\^  -  x(0)>   is  the  circular  correlation  of   the 
((rU)P))      J 

sequence  J  exp  -  j  (2rf/N)r^n^p  \   , 


In  matrix  notation,  Eg. i  1-36]  iDay  be  written  as 


i  i  i 

X   -  x  (0)  D  =  F  x  (n) 


(1-39) 


where 


D  = 


(1-^0) 
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and  X   aid  x  (d)  are  column  vectors  of   size   (N-1)   derived 
frcm  X(k)  and  x  (n)  by  deleting,  respectively,  X (0)  and  x (0)  . 

i 

Ihe   matrix  F    can   be   factored  into    three 

matrices 


•    t 

F   =PCI  (1-<M) 

where    E    is   ac    (N-1)   permutaticn  matrix,   C   is   an 

[  (N-1)  X  (N-1)  ]  circulant  matrix,  and  P   is  the   transpose   of 


lhe  elements  of  the  C  matrix  are 


C     =  F  ,  n,k=1,...,N-1         (1-42) 

r,k      (n)  p,  (k) p 


where  "r"  is  a  primitive  root  of  N  and  "p"  implies  permuted, 
The  elements  cf  the  permutaticn  matrix, P,  are 


E    =  J  n,k=1,...,N-1  (1-43) 

o,k     (n)p,k 


where    C         is   the   Kronecher   delta   function.    The 
S(n)E,k 

circulant  matrix, C,  can  be  implemented  with  a   recirculating 
tisrsversal  iilxei  whcse  tap  weights  are  determined  by 


(n)  p 
h (r)  =  W        ,  n=1,...,N-1  (1-44) 
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Ihus  the  prime  transform  algorithm  offers  more 
sinjlif icaticr  than  the  C2T  algorithm  because  in  hardware 
iif lementaticn  it  is  possible  to  eliminate  the  pre-  and 
pest-multipliers.  The  processing  operations  maj  be 
sunaarized  as 

1.  p ernttaticn  of  the  input  data 

2.  circular  correlation 

3.  per  nutation  of  the  output  Fourier  coefficients 
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HSfiESlBE  AKE  SCFTWA£S  IHELEMEHTATICN 


in  fccth  hardware  and  software  implementation  of  theses 
algcrithns,  the  goal  is  to  have  the  ability  to  do  real-time 
signal  piccessing  for  the  frequency  range  cf  interest.  In 
spectral  acclysis  applications,  signal  prosessing  takes 
place  in  real  time  when  the  spectrum  level  cf  each  frequency 
ccmponent  is  updated  at  a  rate  equal  tc  at  least  the 
bandwidth  p  (Bz)  cf  the  system. 

Depending  upon  the  system,  the  signal  for  processing  may 
he  in  any  cf  three  forms,  i.e., 


1.   Analog-*  toth  time  and  signal   amplitude   are 
ccttinuous 


2.   Sanrled  analog-*  time  is   discrete  and  signal 
amplitude  is  continuous. 


3.   Eicital-*  time  and  signal  amplitude  are   fccth 
discrete 


1-   ii3Jica  Spectral  An^IlJi^ 


Scanning  Analyzers 


The  conventional  scanning-filter  analyzer 
consists  cf  a  single  bandpass  filter,  whose  center  frequency 
is  "shifted"  by  using  frequency- translation  techniques. 
Typically  a  single  bandpass  filter  is  employed  and  the 
filter  analyzes  cne  freguency  at  a  time. 
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Since  a  single  filter  is  employed  and  is  applied 
successively  tc  each  frequency,  it  must  dwell  for  a  period 
of  1/g  at  each  frequency.  Thus  (1/f  )  X  No.  of  steps  equals 
tie  total  analyzing  time  "T".  The  narrower  the  filter  the 
longer  is  the  total  sweep  time.  Since  each  freguency 
component  is  updated  only  once  per  time  Tr  real  time  signal 
prccessirg  is  ret  possible. 

t.   tultifilter  Spectrum  Analyzers 

Cccstact  bandwidth  multifilter  analyzers  were 
developed  for  real  time  analysis.  In  this  type  of  analyzer, 
the  sicnal  passes  through  many  parallel  filters 
simultaneously  ,  typically  500  contiguous  magnetostrictive 
filters  per  analysis  range.  The  filters  and  the  associated 
circuits  must  possess  uniform  and  constant  gain  so  that  a 
fiied  anplitude  sine  wave  of  any  freguency  in  the  analysis 
range  results  in  a  spectrum  output  which  is  constant, 
hchever,  frequency  flatness  of  better  than  ±3dB  is  difficult 
tc  maintain,  thus  the  magnetcstrictive  multifilter  analyzers 
are  net  used  eitensively. 

2»   li^i^Ji  ii?S.ctral  Analysis 

The  all  digital  spectrum  analyzers  employ  a  computer 
or  computer-like  circuitry  which  is  programmed  to  implement 
the  FrT  algorithm.  Some  of  the  structural  factors  fcr  the 
design  cf  such  a  system  are  discussed  helow. 
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a.   Technology 

Easic  speed  and  functional  blocks  available  are 
dictated  ty  technology.  Table-Ill  shows  a  comparison  of  the 
najcr  logic  families.  Because  of  speed-complexity 
tradeoffs,  tie  choice  of  logic  family  is  usually  not  a 
clear-cut  decesion.  This  complexity  encompasses  many  areas, 
tut  important  are  available  MSI  and  LSI  functional  modules, 
system  inter ccnnections,  multiple  function  modules,  and  type 
of  multiplier.  To  gain  speed  in  the  fft  algorithm  a  fast 
multiplier  is  required  in  the  butterflies,  thus  the  type  of 
multiplier  is  an  important  consideration.  Table-IV  is  a 
comparison  of  multipliers. 

TABLE-IV 

MULTIPLIER  COMPARISON   TABLE 


Configuration 

Type  of 
Basic  Package 

Component 
Technology 

Comp 

16  x  12 

nation  1 
(nsec) 
16  x  8 

"ime 
9  :<  9 

Pack 
16  x  12 

age  Cow 
16  x  8 

U 
9x  9 

Fast  trapezoid 

Custom  two- 

ECL 

32 

28 

22 

112 

82 

45 

(Fig.  8.35) 
Tree 

bit  adder 
package 
Commercial 
four-bit 
adder 

ECL 

62 

48 

40 

160 

120 

85 

Tree 

package 

Commercial 

four-bit 

adder 

TTL 

155 

100 

85 

160 

120 

85 

package 
(2  x  4)-bit 

TTL 

261 

195 

135 

80 

58 

30 

(commercial) 

array  pack- 

Two bits  at  a 

age 
Custom  two- 

ECL 

130 

120 

100 

40 

36 

25 

time  add- 

bit  adder 

shift 
Two  bits  at  a 

package 
Commercial 

ECL 

210 

140 

116 

40 

36 

25 

time  add- 

two-bit 

shift 

adder 

Two  bits  at  a 

package 

Commercial 

TTL 

450 

300 

250 

40 

36 

25 

time  add- 

two-bit 

shift 

adder 
package 

Taken  from  Ref.[l7]  . 
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t.   algorithm 

lhe  structure  of  the  FFT  algorithm  is  a  major 
factor  in  determining  whether  or  not  special-purpose 
hardware  cr  a  general-purpose  computer  is  used  in  a 
particular  application.  Table-V  lists  the  arithmetic 
operations  required  for  different  radix  algorithms.  When 
inplemented  with  software,  the  radix-8  algorithm  is  near 
cptimum.  Hc*ever,  it  lacks  the  flexability  required  of  a 
general-cur pcse  ccnputer  and  thus  is  only  used  in  special 
purpose  applicaticns.  Although  the  radix-2  requires 
additional  ecu putations,  its  simplicity  has  made  it  poplular 
for  use  in  general  purpose  computers  were  real  time  analysis 
is  ret  a  requirement. 

c.   Hardware  Architecture 

fchen  inplemented  with  hardware,  the  radix-4 
algorithm  hes  received  the  most  attention.  Hardware 
iaplementaticc  cf  the  FFT  algorithm  has  gained  speed  by 
trading  eff  flexability.  The  speed  is  gained  by  a  cemplex 
nix  cf  several  types  of  parallism.  One  such  type  is 
pipelining,  i.e.,  the  process  is  broken  up'  into  several 
sequential  tasks  and  execution  proceeds  assembly-line  style. 
Ey  pipelining  many  processing  steps  can  be  progress  at  any 
given  tine  and  thus  performing  real  time  analysis.  Given 
that  an  N-pcint  FF1  requires  (N/2) log  N   butterflies,   there 

are  four  hasic  ways  to  organize  the  processor. 
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1.   Sequential-*  all   (N/2) log  N  butterflies  are 
peifcrmed  sequentially. 


2.  Cascade-*   if  log^N  arithmetic  units  (AD)  are 

available,  the  algotihm  can  be   pipelined   sc 

that   each  AU   computes   (N/2)   butterflies 
before  advancing  data  through  the  pipe. 

3.  Earallel-Iterative-*  if  (N/2)  AU»s  are 
available,  the  algorithm  can  be  prccess  in  a 
lcc^N   stage   pipeline.    This  reguires  leg  K 

^  2 

tutterfly  times. 


4.   Array-*  if  (N/2) log  N  AU's  are  available,  the 

entire  array  can  oipelined  into  one  execution 
tine  . 


Table-VI  ccupares  these  different  organizations  fcr  an 
N=1G24  FF1  ate  an  execution  time  of  1  f s  The  total  execution 
time  does  net  take  into  account  the  I/O  overhead  time. 

TABLE-VI 
Comparison  of  Processor  Organizations 


Machine      Arith.  Units      Total 
Organization   (Butterflies)   Execution 

Time  (us) 


Sequential 

1. 

5120. 

Cascade 

10. 

512 . 

Parallel 
Iterative 

512 . 

10. 

Array 

5120. 

1. 

6* 


In  addticn  tc  pipelining,  among  the  types  cf  parallism  are 

1.  Tine  overlap  cf  control  and  memory  functions. 

2.  The   addition  of  a  small  amount  of  high  speed 
meiicry. 

3.  Higher  radix  algorithms. 

Tahle-VII  ccnpares   some   important   parameters   of   a   few 
modern  sicca]  prccessing  systems. 

3«   -ambled  Analog  Spectral  Analysis 

like  the  FFT  algorithm  the  CZT  algorithm  can  ccmpute 
the  freguency  ccmpcnents  of  a  signal  mete  efficiently  than 
by  direct  evaluation  of  the  DFT.  However,  for  relatively 
large  values  of  N,  the  CZT  algorithn  is  less  efficient  when 
cempared  tc  the  FFT  algorithm.  For  this  reason  interest  tc 
digitally  inpleaent  the  CZT  algorithm  waned  while  great 
strides  here  being  taken  to  develcpe  FFT  software  and 
haicware. 

The  CZT  does  offer  many  signal  prccessing 
advantages,  bet  only  if  it  can  be  implemented  using  sampled 
analog  techcigues.  With  the  developement  of  Surface 
Accustic  have  (SAft)  and  Charge  Transfer  Devices  (CTD) , 
interest  in  the  CZT  has  been  re-kindled.  The  prime 
advadvantace  cf  these  devices  ever  digtal  eguipment  is  that 
when  the  ccrtinucus  analog  signal  is  sampled  the  analog 
value  is  retained,  while  in  a  digital  system  the  analog 
value  must  be  quantized  which  is  a  source  cf  error. 

CTE's,  which  include  the  Charge  Coupled  Device  (CCD) 
and   Eucket-Erigade   Device   (SBD) ,   sample   the   optical  or 


62 


electrical  input  and  convert  it  to  a  packet  of  charge  whose 
magnitude  is  prctortinal  to  the  input.  The  charge  packet 
neves  dean  the  CTD  stages,  shift  register  fashion,  under  the 
ccntrcl  cf  clock  voltages,  The  charge  packets  are 
transferred  frcn  one  stage  to  the  next  until  they  readh  an 
output  stace  where  the  magnitude  of  the  charge  packet  is 
sensed  anc  converted  to  an  electrical  anlacg  signal.  The 
clcck  frequency  and  the  number  of  stages  transferred 
determine  tte  time  delay  introduced. 

In  c  SAW  device  a  mechanical  ripple,  or  wave, 
travels  alone  a  solid  polished  surface  in  much  the  same 
manner  as  an  ocean  wave  travels  through  water.  Aluninum 
electrodes  are  usually  deposited  on  a  piezoelectric  material 
such  as  £1  quartz,  and  they  act  to  launch  and  receive  the 
accustic  fcaves.  These  transducers  convert  energy  between 
electrical  ard  mecharical  domains. 

In  a  CCD  the  delay  line  is  tapped  and  weighted  by 
using  a  split-gate  weighting  technigue  which  is  incorporated 
intc  the  mask  at  the  time  of  fabrication  or  ty  weighting  the 
taps  with  external  resistors.  The  tap  weight  in  a  SAW 
device  is  photc-lithographically  determined  and  is 
prcpcrticral  to  the  length  of  the  interdigitaticn  which 
accepts  the  accustic  wave.  The  pre-and  post-multiplies  are 
accomplished  using  wideband  four  analog  multipliers  with 
weights  generated  ty  programmable  read  orlj  memories  (EECM) 
and  ccucled  through  D/A's  or  by  multiplying  D/A's  with  EROM 
weight  stcrace. 

Ihe  transversal  filter  implementation  of  the  CZT 
algcrithi  achieves  its  high  speed  not  by  reducing  the  cumber 
cf  multiplies  like  the  FFT,  but  by  putting  the  eguivalent  of 
M  parallel  multipliers  all  on  the  same  chip.  This  reduces 
the  number  ci    computational  steps  from  Nlog^N  for  the  IM  to 
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N  fcr  the  CZ1.  Also  the  N  additons  per  component  are 
sinultanecus.  Currently,  the  CCD  CZT  can  operate  in  real 
tine  at  speeds  up  to  5MHz  and  the  EflD  CZT  is  limited  to  a 
sample  rate  cf  a  feu  hundred  KHz.  For  higher  sample  rates, 
up  to  several  fifiz,  SAW  devices  are  attractive.  In  addition, 
the  sampled  analog  devices  require  lew  power  and  are  light 
weight.  The  advantages  of  the  digital  FFT  are  its 
flexability  and  accuracy.  Programability  is  teing  developed 
for  the  CCD  and  SAW  filters  with  up  tc  128  digitally 
switched  taps  have  been  built. 

As  rcted  earlier,  the  prime  transform  effers 
additional  simplicity  by  elimininating  the  multipliers 
reguired  in  the  CZi  implementation.  The  multipliers  are 
replaced  with  analog  permuter  memories.  The  commercially 
available  serial  access  memory  stores  analcg  samples  as 
charges  in  an  array  of  MOS  capacitors  under  the  ccntrcl  of 
read  in  and  read  out  shift  registers.  Currently  the  maximum 
sample  rate  fcr  this  type  of  device  is  5  MBz. 
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12-  IAM1   IOUBIJB  TBANSFOBM  ALGCglTHMS 


A.   MATH! EPICAL  ILLUSTRATION  OF  THE  DECIMATICN-IN-TIME 
(EII)  ALGCfllEfi 


lo   illustrate   the   FFT  algorithm,  the  number  cf  sample 

y 
pcints  of  x(n)  is  assumed  to  te  a  power  of  2,  i.e.,  N  =   2  . 

In  this  exanple  N=4=22  andV=2. 


For  K=4,  Eg. [1-3]  can  be  written  a; 


X(0)  =  x(0)W°  +  x(l)W°  +  x(2)W°  +  x(3)W° 
X(l)  =  x(0)W°  +  xiDW1  +  x(2)W2  +  x(3)W3 


X(2)  =  x(0)W°  +  x(l)W2 


x(2)W^  +  x(3)W6 


X(3)  =  x(0)W°  +  x(l)W3  +  x(2)W6  +  x(3)W9 
or  written  in  matrix  form 


(2-1) 


x(o) 

X(l) 

X(2) 

X(3) 

w°  w1 


o   3 
w   w 


w   w 


w   w 


o" 

"x(oj" 

3 

x(l) 

6 

x(2) 

9 

x(3) 

(2-2) 
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The  first  step  in  developing  the  FFT  algorithm  for  this 
example  is  tc  rewrite  Eg. [2-2]  using  the  following  relation 
fica  number  tfcecrj[23]. 


e)<     rk  mcd/N) 


(2-3) 


where  [nk  mcc(N)  ]  is  the  remainder  upon  division  of  nk  by  N, 
Hence  if  1=4,  n=3,k=3 


S     1 

e  =  h 


(2-4) 


since 


nk     S 
i    =  I   =  exp[  (-j2n/4)  (9)  ]  =  exp(-j9rT/2) 


=  €xp(-jTT/2)  =  exp[  (-J2T/4)  (1)  ] 


1     rk  mod  (N) 

=  w  =  i 


(2-5) 


Thus  the  new  uatrix  is  written  a; 


X(0) 
X(l) 
X(2) 
X(3) 


L1 


i 


1      1 

~x(0)" 

w2    w^ 

x(l) 

w°    w2 

x(2) 

2         1 

vr    vr 

x(3) 

(2-6) 


The  second  step  is  the  factorization  of  the  sguare 
matrix  in  Ig.[2-€].  As  an  aid  in  showing  hew  Eg. [2-6]  is 
factored,  a  new  notation  is  introduced  at  this  time,  where 
the   integers  o  and  k  are  now  represented  as  binary  nutrbers. 
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n    =    0,1,2,3  — >   n  (n    ,n    )     =    CO,    01,    10,    11 


k    =    C,1,2,3   -->   k(k    ,k   )     =    GO,    01,     10,     11 
1      0 


Then    writing    the    case    10   equivalents 


n=2n      +n  ,k=2k      +    k 

10  10 


(2-7) 


Eg. [1-9]   can   new   fc€    written   as 

1  1 

X(k. 


(2n,+nn)(2kn+kn) 


V^      V  C2n1+nnH2k1+kn 

:1^0)    =  ^x(ni,n0)W        10  10 


(2-8) 


nQ=0  n,=0 


Tt€  important  pemt  here,  is  that  the  single  summation  of 
Eg. [1-9]  mest  te  replaced  by  "  "  summations  in  order  to 
enumerate  all  the  tits  of  the  binary  representation  of  n. 

Continuing 


w(2n1+n0)(2k1+k0)  =   (21^+^)2^  w(2k1+k0)n0 


'  ^n  k' 

w  i  i 


w2Vi  w(2ki+Vno 


w2koni  w(2ki+ko)no 


(2-9) 


The  tern  in  tiackets  is  equal  to  unity  sines 


A"lkl 


niki 


j2  V+' 


nlkl 


nlkl 
1  L   L   =  1 


.2-10 
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Eg. [2-8]      cat    now    te    written    as 


X( 


± 


nQ=0 


£ 


x(n1,nQ)W 


2k„n 


0^1 


Uk1+k0)n0 


(2-11) 


Each   of   the   summations   are  treated  individually  and  lead 
directly   tc   the   correct   factorization   cf   Eg. [2-6] 
Eehritinc  the  bracketed  summation 


x1(k0,n0)  =  )  x(n1(n0)  W 


2Vl 


n1=0 


Directly  evaluating  this  summation  gives 


(2-12) 


x1(0,0)  =  x(0,0)  +■  x(l,0)W( 
x1(0fl)  =  x(0,l)  +  x(l,l)WC 
x1(l,0)  =  x(0,0)  +  x(l,0)W 
x1(l,l)  =  x(0fl)  +  x(l,l)W; 
re-writing  in  matrix  notation 


(2-13) 


x1(0,0) 

1 

0 

w° 

0 

~x(0,0)" 

x1(0,l) 

0 

1 

0 

w° 

x(0,l) 

x1(l,0) 

1 

0 

w2 

0 

x(l,0) 

x1(l,l)_ 

0 

1 

0 

w2 

_x(l,D_ 

(2-14) 
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Now  writing  the  outer  summation  as 


x2( 


Xl(k0,nQ)  W 


(2k1+k0)n0 


nQ=0 


(2-15) 


Directly  evaluating  the  summation  gives 


x2(0,0)  =  x1(0,0)  +  x1(0,l)W° 

x2(0,l)  =  x1(0,0)  +  x1(0,l)W2 

x2(l,0)  =  x1(l,0)  +  x1(l,l)W1 

x2(l,l)  =  x1(l,0)  +  x1(l,l)W3 


(2-16) 


Then  in  matrix  notation 


x2(0,07 

1 

w° 

0 

0 

^(0,0)" 

x2(0,l) 

1 

w2 

0 

0 

x]_(0,l) 

x2(l,0) 

— 

0 

0 

1 

w1 

x1(l,0) 

x2(l,l) 

0 

0 

1 

w3 

x1(l,l) 

(2-17) 


Ccntininc    Eg£.[2-12]    and   [2-15]    results    in 


'S'V    =    W    V 


(2-18) 


or  hrititg  ir  latrix  form 


x(o) 

"l 

W°   0 

0~ 

l 

0 

w°  o" 

'x(oT| 

X(2) 

1 

2 

vr  o 

0 

0 

1 

0   w° 

x(l) 

X(l) 

0 

0    1 

w1 

1 

0 

W2   0 

x(2) 

X(3) 

_0 

0    1 

w3 

_0 

1 

2 
0    VT 

_x(3l 

(2-19 
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A   signal    flew    graph    of   Eg. [2-19]    is    shewn    in   Fig. [13], 


COMPUTATION  ARRAYS 


x2(0) 


x(o) 


*2m   X(2) 


*--^x2(2)   x(l) 


X(3) 


Jigure  13  -   FFT  Signal  Flcwgraph,  N=4 . 


Note  that  the  final  results  are  in  bit  reversed  order   with 

respect  tc  tfce  disered  values  X(k  ,  k  ).   This  is  simply  the 

10 

scrambling   which   results   frcm   the   FFT   algorithm.    The 

results  are  easily  unscrambled  by  bit  reversing,  i.e., 


X(CC)  X  (00) 

X(1C)  flips  or  reverses  to    X(01) 

X(C1)  X(10) 

X(11)  X(11) 


(2-20) 


The  algorithm  just  described  is  known  as  the 
Cccley-Tukey  Agcrithm  or  decimation-in-time  (EIT)  algorithm, 
since  at  each  stage  cf  the  process  the  input  sequence  (i.e., 
time  secuence)  is  divided  into  smaller  sequences  for 
processing . 

The  hasic  building  block  of  the  algorithm  is  the 
"fcutterflj"  as  shefcn  in  Fig. [14]. 
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AV         PX   =  A+ 

es  \y=a-bi 


Figure  14  -   Decimation-In-Time  Eutterfly. 

Each  butterfly  requires  one  complex  multiply  (four  real 
multipliers)  and  two  complex  adds.  The  inportance  cf  the 
butterfly  structure  to  the  FFT  flew  graph  is  that  only  one 
additional  nemory  lecatien  is  required  tc  transform  an 
N-pcint  seguence  stored  in  memory.  Thus  the  intermediate 
results  cf  the  FFT  are  stored  in  the  same  locations  in  which 
the  original  data  being  transformed  was  stored.  Since  this 
algcritha  stores  both  the  input  and  output  sequences  in  the 
same  location,  it  is  called  an  in-place  algorithm. 

As  pointed  out  earlier  the  results  of  the  DIT  algorithm 
are  bit  reversed.  This  means  that  to  get  a  normal  ordered 
output,  the  input  data  must  first  be  shuffled.  Figs. [15a] 
arc  [15b]  present  the  canonic  signal  flew  graphs  for  the  DIT 
algorithm  for  the  case  N=8.  Note  that  the  data  is  naturally 
ordered  ir  Fig. [15a]  and  is  in  inverse  order  in  Fig. [15b]. 
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B.   EFCIE21ICK-IN-fR£QUENCY  (DIP)  ALGORITHM 


Anctter  popular  FFT  algorithm  is  the  Sande-Tukey 
Algorithm  or  decimation-in-f reguency  (EIF)  algorithm.  It  is 
very  siuilai  to  the  DIT  algorithm  but  differs  in  the 
following  wa5s. 


1 . 


lhe 
two 
tie 


input  seguence  fx(n)}  is  partitioned  intc 
twc  sequences  each  of  length  (N/21  samples 


following   manner 


The   first 


IE 

seguenc 


Jx  (n)}  consists  of  the  first  (N/2)  points  of 
|x(e)},  where  as  the  second  seguence  {x  (n) } 
consists  of  the  last  (N/2)  points  of  {x(n)}. 


For  a  given  case  in  the  DIT  algorithm  the 
incut  is  tit-reversed  while  the  output  is  ic 
natural  order,  whereas  the  reverse  is  true  for 
the  EIF  algorithm.  In  general,  however.  hctfc 
the  EIT  and  EIF  can  go  from  normal  to  shuffled 
data  and  vice  versa. 


3.  lie  EIF  butterfly  is  slightly  different  as 
shc*E  in  Pig.[l6j.  Ihe  difference  fceing  that 
the   complex   multiplication  takes  place  after 


the   add-suttract 
algorithm 


operation 


the 


Elf 


figure  16  -   Decimaticn-In-Freguency  Eutterfly 
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Taking  a  lock  at  the  similarities,  tcth  algorithms 
require  ct  the  crd€r  of  (N  log  N)  operations  tc  compute   the 

DPI.    Ectb  algorithms  can  be  done  in-place  and  both  need  tc 

perfcrm  tit  reversal  at  some  place  during  the  computation. 

Figs. [15c]  and  [15d]  present  the  cancnic  signal  flow 
craphs  fcr  the  DIF  algorithm  fcr  the  case  N=8.  The  data  is 
natcrallj  ordered  in  Fig.£15c]  and  is  in  inverse  crder  in 
Fig.«15d].  Ihe  two  most  effective  metheds  are  those 
illustrated  in  Figs. [15b]  and  [15c]  since  they  provide  the 
peters  cf  U  ir  the  correct  crder  needed  fcr  computation. 
This  eliminates  the  need  for  storage  tables. 

The  III  algorithms  examined  above  by  nc  means  eahaust 
all  the  possibilities.  Numerous  algorithms  have  been 
developed,  all  kith  cne  basic  purpose  in  mind,  and  that  is 
tc  reduce  the  computation  time. 
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EXAMEIES  CSING  THE  FFT  ALGORITHM 


As  a  test  cf  the  EFT  algorithm  to  compute  the  DFT,  CASES 
1,2,  and  3,  as  described  in  Appendix  A,  were  used  as  inputs 
for  N  =  32.  For  the  EIT  algorithm  used,  N  is  reguired  to  be 
highly  ccmpcsite  and  a  power  cf  2.  The  value  of  32  was  used 
so  that  the  results  of  this  algorithm  could  te  compared  to 
the  results  cf  the  CZT  and  prime  transform  algorithms.  In 
the  latter  twc  algorithms  N=31.  All  figures  are  drawn  with 
normalized  frequency,  but  for  the  sake  of  discussion,  assume 
that  the  incut  sinusoids  have  freguency  in  Hertz  anc  the 
rectangular  cata  window  "T"  is  egual  to  1  second. 

The  fundamental  freguency  is 


E   =1/1=1  BZ 

c 


which  is  the  ninimum   recognizable   freguency   component   or 
resolution,   the  analyzing  bandwidth  is 


?     =  (K/2)  (I  )  =  16  Hz 

max  o 


The  output  is  a  one-sided  spectra  consisting  of  [|N/2)  +  1] 
freguency  components.  Figs. [17]  and  [18]  are  the  results 
for  CASE  1  and  those  for  CASE  2  are  shown  in  Fig.[1S]  and 
[2C].  In  the  results  cf  CASE  3,  the  22  Hz  component  aliases 
as  a  10  E2  component. 
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III.   CHIBP-Z-TEANSJCRM  COMPUTER  SIMULATION 


A.   CZT  IFOGfAIT  DEVEICPEMENT 


A  computer  program  was  written  to  sinulate  the  CZT 
algcrithc  fcr  computing  the  DFT  of  an  arbitrary  waveform. 
This  waveform  nay  te  real  or  complex.  In  carrying  cut  the 
piccessinc  operations,  the  real  and  imaginary  components  are 
represented  fcj  subscripts  "r"  and  "i"  respectively.  The 
simulation  program  is  a  result  of  the  following  derivation. 
Proceeding,  let 


then 


x-^  input  sequence 


-n  n2/2  „  „ 

f  =  A   W    — >  complex  pre-multiplier         (3-1) 


y  =  >f  (3-2) 


y  =  x  f   ♦  jx  f   +  jx  f   +  j2x,f.  (3-3) 

rrri     ir      11 


simplif yirg 


y   =  x  f   -  x  f  (3-4) 

r    r  r    i  i 


y   =i(xf   ♦  x  f  )  (3-5) 

i      r  i    i  r 
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Now   letting 


v   =    fi 


n2/2 


g    =    j    *    v 


(3-6) 
(3-7) 


where    •■*"    denotes   convolution 


g    =    y    *v      ♦   y    *jv      +    jy    *v      +    jy    *jv 
rr  ri  ir  ii 


(3-8) 


g      =    y   *v      -   y   *v 

r  r     r  i     i 


(3-9) 


g        =    y    *jy        +     jy     *V 

i  r        i  i      r 


(3-10) 


Id   the   program   this   convolution   is    simulated    as   a 

transversal   filter  with  2N-1  conplex  taps  v       to  v 

-<N-1)       (N-1) 


-ez/2 
v   =  H         ,  n  =  -  (N-1)  to  (N-1) 

n 


(3-11) 


ana 


i  =  €Xf  (-J2TT/E) 

The  third  uajcr  step  in  the   process,   that   of   complex 
pcst-multip licatioc,  is  performed  fcy  letting 


d  =  K 


W 


(3-12) 
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X  =  gd 


X=gd   +  jg  d   +  jg  d   +  j2g  d 
r  r     ri     ir      ii 


(3-13) 
(3-14) 


simplifyicg 


X   =  g  d   -  g  a 
r    r  r    i  i 


(3-15) 


x.  s  j(g  d,  +g.d) 
i      r  i    i  r 


(3-16) 


Fcr  the  special  case  of  the  DPT,  A  =1  and  9  =0  in  Eg. [1-23], 

o        c 

~n 
«  =1   and  0   =-(1/N)   in  Eg. [1-24].   Thus  A   =1,  W  =  e*p  (-  j2 
c  c 

tt/N)  and  Ec.[3-1],  [3-6],  and  [3-12]  can   be   thought   cf   as 

complex    eifcnential   seguences   of   linearly   increasing 
frequency.   They  are  represented  as  follows: 


f   =  cos  (TTn2/N) 

r 


(3-17) 


f   =  ^sir  (TTn2/N) 
i 


(3-18) 


V    =  COS  (TTn2/N) 


(3-19) 


v   =  sin  (rrn2/N) 
i 


(3-20) 


d    =  CCS  (TTk2/M) 

r 


(3-21) 


d   =  -sin  (TTk2/M) 

i 


(3-22) 
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Substituting  Eg. £3-17]  thru  [3-22]  into  Egs.[3-4],  [3-5], 
[3-9],  [3-1C:#  [3-15],  and  [3-16],  respectively,  describes 
the  implementation  of  the  CZT  as  shown  pictorially  in 
Fig. [23]. 


I      =   x   ccs(TTn2/N)    +   x   sin(TTn2/N) 
r  c  i 


(3-23) 


y.    =   j[-x   sin(TTn2/N)    +    x  cos  frrnZ/N)  ] 


(3-24) 


c      =   [y    *ccs(TTn2/N)  ]    -   [y    *sin  (TTn2/N)  ]  (3-25) 

r  r  i 


9.    =   [y    *jsin  (TTn2/N)  ]   +    [jy    *cos  (TTn2/N)  ]  (3-26) 

ii  i 


X      =   g   ccs.f^/M)     +    g   sin(TTk2/M) 
r  r  i 


(3-27) 


X      =    j[-g    sin  (tt!c2/B)     ♦    g    cos  (TTk2/M)  ] 

i  r  i 


(3-28) 


lh€      simulation      program      is      listed      in   Appendix   D    anc    is   a 
working    icdel    of   Fig. [23]. 
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EXAflEIES  USING  THE  CZT  ALGORITHM 


As  a  test  cf  the  CZT  algorithm  to  compute  the  DPT,  CASES 
1,2,  and  3  as  described  in  Appendix  A,  were  used  as  inputs 
for  N=31.  Ike  imaginary  component  of  the  input  is  set  equal 
to  zero  since  the  input  consists  of  only  real  data.  All 
figures  are  crawn  «ith  normalized  frequency,  but  for  the 
sake  cf  discussion,  assume  that  the  input  sinusoids  have 
frequency  in  Bertz  and  the  data  window  "T"  is  equal  to  1 
second.  The  data  window  used  for  all  cases  is  rectanqular. 
Thcs  the  fuccamental  frequency  is 


1   = 


=  1  Hz 


(1)    N(T  )     (31)  (1/31) 

s 


lihich  is  ic  effect  the  minunum  recognizable  frequency 
ccjrpcnent  ci  resolution.  Ihe  analyzing  bandwidth  or 
frequency  racce  is 


E     =  (N/2)  (E  )  =  15.5  Hz 

iax  o 


However,  since  the  resolution  is  1  Hz,  the  output  components 
consist  cf  the  D-C  component,  the  fundamental,  and  the  first 
14  harmonics  for  a  frequency  range  of  15  Hz. 

As  ic  the  EFT  algorithm,  the  output  of  the  CZT  algorithm 
is  a  one-sided  specta  consisting  cf  [  (N/2) +1  ]  integer 
frequency  ccapcnents.  Figs. [24]  and  [25]  are  the  expected 
results  for  CASE  1  and  Figs. [26]  and  [27]  are  the  results 
for  CASE  2.  In  the  results  of  CASE  3,  it  should  he  noted 
that   the   aliasing   frequency   of   22   Hz  folds  back  as  a  9 
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Hz, whereas  in  the  FFT  it  folded  tack  as  a  10  Hz  ccmpcnent. 
This  occurs  tecause  N  is  odd  in  the  CZT  where  as  in  the  FFT, 
N  is  required  to  be  even  and  highly  composite.  For  the  FFT 
alccritha  N  =22. 
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IV.  IMM    TMNSFORM  COMPUTER  SIMULATION 


ERIME  TRANSFORM  COMPUTER  EBOGRAM  DEVELCEEMENT 


So  as  tc  nake  the  mathematics  cleat,  a  simple  example  is 
given   tc  illustrate  the  required  matrices  cf  the  algorithm. 

Consider  the  case  N  =  5  and  let  9  =  exp(-j2  tt  /N)  . 
Evaluating  the  DFT  equation  directly  results  ir.  the 
following  N  eguaticns. 


X(0)  =  x(0)W°  +  x(l)W°  +  x(2)W°  +  x(3)W°  +  x(^)W° 
X(l)  =  x(0)W°  +  xCDW1  +  x(2)W2  +  x(3)W3  +  x(^)W^ 
X(2)  =  x(0)W°  +  x(l)W2  +  x(2)W^  +  x(3)W6  +  x(^)W8   (>-D 
X(3)  =  x(0)W°  +  x(l)W3  +  x(2)W6  +  x(3)W9  +  x(^)W12 
X(k)    =  x(0)W°  +  x(l)W^  +  x(2)W8  +  x(3)W12+  x(^)W16 
Written  in  matrix  form 

x(k-) 


Fj  Lx(n)j 


X(0) 
X(l) 
X(2) 
X(3) 
XWJ 


wu 
w3 
w6 


w 


,12 


o" 

x(0) 

k 

x(l) 

r8 

x(2) 

^ 

x(3) 

16 

_x(^)_ 

(^-2) 
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The  matrix  {F}  is  written  using  the  relationships 


rnk 


W1^  mod  N 


W°  =  1 


For  example  w9  =  w^  since  for  n=3,k=3 

W«*-¥*-  exp[(-^)(9)]  =exp(-jl|* 
=  exp(-j|2T)  =  exp[(-j^)(Z|)]  =  w4 

therefore 


(4-3) 


x(o) 

ri 

i 

1 

1 

X(l) 

i 

w1 

w2 

w3 

X(2) 

= 

i 

w2 

w^ 

w1 

x(3) 

i 

W3 

w1 

vA 

x(4) 

i 

w* 

w3 

w2 

W" 


x(0) 
x(l) 
x(2) 
x(3) 


(4-4) 


Then  deleting  X(0)  and  x(0)  ,  Eq.  4-4   is 


"x(i)" 

1 

w1 

w2 

w3 

w^ 

x(l) 

X(2) 
X(3) 

-  x(0) 

1 

= 

w2 
W3 

w1 

w1 

w3 
w2 

x(2) 
x(3) 

_X(4) 

A 

w^ 

W3 

w2 

w1] 

x(4) 

■ 

X 

D 

t 

p 

x(n) 

(4-5) 


For  N=5  and  r=2 ,  there  is  a  one  to  one  mapping  of 
integers  n,k  to  the  integers  (n)p,(k)p,  i.e., 

21  =  2  mod  5,  22  =  k   mod  5,  23  =  3  mod  5,  2^  =  1  mod  5 


the 


n,k 

1 

2 

3 

4 

(n)p,(k)p 

2 

4 

3 

1 

(4-6) 
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Matrix  [F  ]  is  factored  into  three  matrices. 


F  =  p^CP 

The  elements  of  the  [C]  matrix  are  found  by  using 
relationship  given  by  Eq.fl-tel.  Thus 


the 


C  = 


w3 


WJ 


,2    ,.M 


w1 

w2" 

w2 

w^ 

w^ 

W3 

w3 

w1 

W"   W 

The  elements  of  the  permutation  matrix  are  given 
Eq.fl-^l. 


(4-7) 


by 


0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

0 

0 

0 

(4-8) 


The  elements  of  the  transpose  of  P  are 

0    1 

0  0 

1  0 
0    0 


0 

0 

1 

0 

0 

0 

0 

1 

(4-9) 


A  functional  diagram  of  the  prime  transform  algorithm 
for  computing  the  Fourier  coefficients  of  real  data  is  as 
shown  in  Fig. [30  J.  The  switch  is  in  the  up  position  for  the 
first  data  sample  and  is  down  for  the  remaining  (N-l)  sam- 
ples. 
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EXAHEIBS  GSING  THE  PRIME  TBANSFORM  ALGORITHM 


Tc  simulate  the  operations  of  the  prime  transform 
algcrithn,  a  ccmputer  program  was  written  for  the  odd  prime, 
N=31.   A  listing  of  the  program  is  given  in  Appendix  E. 

Once  again  Cases  1,2,  and  3  were  employed  in  order  tc 
evaluate  the  prime  transform  algorithm.  The  fundanental 
freguency  and  freguency  range  are  the  same  as  that  cf  the 
CZT  algcritha  for  H=31.  The  results  are  shewn  in  Figs. [31] 
thru  [36].  fls  can  be  seen,  the  results  are  the  same  as 
these  cciputed  by  the  CZT  algorithm. 
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V.   CZT  SENSITIVITY  ANALYSIS 


This  chapter  addresses  the  sources  cf  error  ir  toth 
digital  and  sampled  analog  systems  with  special  attention 
given  tc  the  sampled  analog  implementation  of  the  CZT 
algcrithB. 


A.   DIGI11I  IIT  SYSTEM 


In   a   digital   FFT   system,  the  sources  cf  error  nay  b< 
sumnarized  as  fellows; 


1.   Quantization   of   the   data   at  the  input  A/I 
converter. 


2.   Inaccurate   transformation   due   tc  errors  in 

k 

firite  representation  cf  the  coefficients   w 

N 

in  the  butterflies. 


3.  Eccndoff  noise  in  truncating  the  result  of  a 
multiplication  and  scaling  the  data  tc 
prevent  cverflow. 


An   las   quantization  error  in  the  A/E  converter  can  be 
defined  as 


-b1 
C   =  2    /  12  (5-1) 

e 
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where  1:1  is  the  nuiber  of  bits  used  to  represent  the  input 
data.  This  ncise  dees  not  scale  with  signal  size  and 
dcainates  enly  at  low  signal  levlels. 

The   errcr  due   to   the  inexact  representation  cf  the 

k 

multiplier  coefficients  W   can  be  expressed  as  the  ratio   of 

N 

mean  sguare  cutput  error  to  mean  sguare  cutput  signal  to 
give 


(2)  (t3) 


where  b3  is  the  number  of  bits  used  to  represent  the  FFT 
ccef fici€nts .  This  eguation  predicts  that  the  variance  of 
the  errcr  increases  at  a  rather  slew  rate  with  N. 

A  measure  of  the  rms  noise  cutput  to   the   rms   signal 

cutout   due   tc   errors   produced  by  roundoff  and   signal 

truncation  tc  prevent  overflows  in  the  butterflies  is  given 
by 

ris(eircr)      V^2   ^<0-3)  8  <5-3) 


rms  (signal)       rms  (input) 


where  fc2  is  the  nunber  of  bits  used  to  represent  the  result 
cf  the  operation.  This  upper  bound  of  the  errcr-to-signal 
is  seen  tc  increase  as  "\f¥,  and  is  the  mest  important  source 
cf  error  in  a  digital  FFT. 
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SA&BLIl    ANALOG  CZT  SYSTEM 


Instead  cf  errors  resulting  from  registers  containing  a 
finite  wcrd  length,  as  in  the  digital  FFT  system,  the  cutput 
accuracy  cf  the  sampled  analog  CZT  system  is  limited  by  the 
CCI's  and  Soft's  sensitivity  to  the  environent  (temperature, 
humidity,  etc.,)  and  fabrication  imperfections.  The  sources 
cf  error  cf  a  CCE  CZT  system  using  pre-  and  post-multiplying 
coefficients  stored  in  either  a  ROM  or  PROM  are  as  fellows 

1.  Thermal  noise-*  this  source  cf  error  is 
analogous  to  quantization  in  a  digital  FF1 
because  it  generates  an  error  which  is 
independent  of  signal  level. 


2.   lh€   errcrs  introduced  by  inaccuracies  in  the 
pre-and  post-multipliers. 


3.   Ibe   errors  introduced  by  the  inaccuracies  in 
the  tap  weights  of  the  transversal  filter. 


4.   Th€   errcr  due  to  mismatch  and  nonlinearities 
in  the  differential  current  amplifiers. 


5.   The    error    due    to    mismatches   of   the 
pcst-multiplier  outputs. 

The  encrs  addressed  in  this  study  where  those  due  to 
inaccuracies  in  the  pre-  and  post-multipliers  and  the 
inaccuracies  in  the  tap  weights  cf  the  transversal  filter. 
The  errors  in  the  pre-  and  pcst-multipliers  were  modeled  as 
being  a  percentage  "p"  of  the  theoretical  coefficient  value 
m (n)  ,  i.e. , 

V=  p[m(n).]  (5-4) 
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and  to  te  uniformly  distributed  on  (-V2,  V/2)  .  The 
transversal  filter  was  modeled  as  a  tapped  analog  delay  line 
iith  the  tap  weights  determined  by  the  value  of  the  external 
resistors.  The  error  in  the  value  of  the  tap  weight  was 
modeled  as  being  a  percentage  »p»  of  the  theoretical  tap 
eright  value  h  (n) ,  i.e., 

P  =  ECMb)J  (5-5) 

acd  to  be  uniformly  distributed  en  (-£/2,  b/2)  .  These  error 
models  were  included  in  the  simulation  program  [Appendix  D]. 


SENSI1IVITI  TESTS 


Three  different  input  waveforms  were  used  to  evaluate 
the  perfcrmace  of  the  CZT  algorithm  for  the  following 
percentages 

p  =  CI,  43,  10%,  20%,  3051,  and  40* 

In  all  ttree  tests  N=512  and  a  rectangular  data  window  equal 

tc  4  seconds  was  used.   The  value  of  4  seconds  was  used   for 

ease   of   ^letting   the  output  spectrum.   Thus  the  frequency 

resolution  is  f/f   =  0.25  and  the  frequency  range  is  f/f   = 
s  s 

64. C. 

The  first  waveform  consisted  of  the  impulse  fuECtion 
shewn  in  Jig. [37].  The  time  impulse  responses  of  the  ccsine 
and  sine  filters  of  the  convolver  for  the  different 
percentages  are  shown  in  Figs. [38]  thru  [50].  For 
p=C2,Fig.[51 j  shows  the  expected  results,  i.e.,  equal 
magnitude  freguency  components  covering  the  entire  freguency 
ranee.  As    tie   percentage  of  error  increases,  examination  of 
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the   figures   reveals  that  the  average  vale  cf  the  magnitude 
of  the  frequency  components  is  greater  than  1.0. 
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Fig. [57]  was  the  second  waveform  used  fcr  computing  the 
DFX.  The  nagnitude  of  the  one-sided  spectra  is  shewn  in 
Figs. £58]  thru  [63].  Degradation  of  the  output  spectrum  was 
never  sc  severe  that  the  (sin  x) /x  function  was  not 
recognizable,  hut  as  pointed  out  above,  the  average  value  of 
the  magnitude  sho*s  a  moderate  increase  as  the  value  cf  "p" 
increases. 

The  third  input  waveform,  shown  in  Fig. [64],  consisted 
of  the  fcllcfcirg  components. 


Amp. 

f/f 

s 

Phase 
Shift 

sine 

1.C 

14.125 

0o 

ccsire 

-0.6 

15.0 

0° 

sine 

1.C 

16.0 

00 

cosine 

0.6 

16.25 

0° 

sine 

-o.  u 

17.5 

600 

The  normalized  one-sided  spectra  and  phase  plots  are 
shewn  in  Fig. [65]  thru  [76].  In  Fig. [65],  the  spreading  of 
the  14.125  component  is  due  to  leakage,  the  16.0  and  16.25 
components  appear  to  be  indistinguishable,  however,  this  is 
a  Hiss-leading  result  caused  by  the  continuous  curve  plot. 
If  a  point  graph  had  been  used,  the  two  components  would  be 
readily  identified,  an  expected  result  since   the   freguency 

resolution  is  f/f   =  0.25. 

s 

Examination  of  the  phase  plots  reveals  that  the 
degradation  is  mere  severe  than  that  of  the  magnitude  plots. 
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even  for  snail  values  of  "p".  As  in  the  previous  exanples, 
the  magnitude  plots  also  exhitit  a  general  increase  in  the 
component  magnitude  except  for  the  f/f   =   16.0   component. 

This   component   shows  a  very  distinguishing  decrease  as  "p" 

increases.   Ej   re-examining   the   magnitude   plots   of  the 

impulse   waveform,   it  is   obvious  that   the   f/f   =  16.0 

s 

component  will  sho-w  a   decrease  as  "p"   increases.   These 

plots  also  show  that  a  more  severe  degradation  of  the 
one-sided   spectra   would  result   for   a   similar   set   of 

ccipcnents  near  f/f   =  48.0. 
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VI.   CONCLUSIONS  AND  fiECCflHEMEATIONS 


The  fJT  clgcrithm  is  a  well  established  method  for 
confuting  the  DF1.  The  chirp-z-transf crm  algorithm  can  also 
compute  the  CFT  and  since  the  bulk  cf  the  required 
confutations  can  be  performed  by  a  transversal  filter,  it 
lends  itself  naturally  to  implementation  with  sampled  analog 
devices.  The  prime  transform  algorithm  has  demonstrated  its 
ability  to  compute  the  DFT  and  like  the  CZT  algorithm,  can 
te  iaplenerted  with  sampled  analog  devices.  The  fact  that 
the  multipliers  reguired  by  the  CZT  algorithm  are  replaced 
by  permutinc  memories  in  the  prime  transform  algorithm  may 
or  may  net  te  an  advantage.  It  will  depend  upon  the  speed 
and  dynamic  range  required  for  a  particular  application. 

This  paper  evaluated  the  errors  in  the  weights  cf  a 
transversal  filter  for  weighting  provided  externally  to  the 
CCE  ty  resistors.  Preliminary  results  indicate  that  the  tap 
weights  nay  vary  as  much  as  20%  before  the  freguency 
components  are  greatly  distorted.  However,  additional 
aralysis  needs  to  te  conducted  before  conclusive  results  are 
presented.  In  addition  to  resistor  weighting,  a  split-gate 
weight  technigue  has  teen  developed  which  weights  the  taps 
of  the  transversal  filter  during  fabrication.  For  this  type 
cf  weighting,  quantization  errors  are  introduced.  Ibis 
suggests  that  a  different  error  model  will  re  reguired  to 
conduct  a  setsitivity  analysis. 
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APPENDIX  A 
ALG0BI1HM  TES1  DATA 


Several  test  waveforms  were  computer  generated  to 
evaluate  the  performance  of  the  various  algorithms  used  to 
ccnpute  the  Ell. 

Case  1->  This  waveform  consisted  cf  the  following 
ccapcnents: 


Amp  . 

f/f 

s 

Phase 
Shift 

E-C 

0.8 

0.0 

00 

cosire 

1.0 

2.0 

0o 

sice 

2.0 

6.0 

0O 

ccsiEe 

1.0 

14.0 

450 

Tfce  intert  cf  this  case  was  to  see  how  well  the  algorithms 
identified  the  uagnitude  and  phase  information  cf  each 
freguencv  ccnpcnent. 
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Case-*   This   waveform   consisted   of    the    following 
ccnpcnents 


Amp. 

f/f 
s 

Phase 
Shift 

E-C 

c.e 

0.0 

0o 

cosine 

1-C 

2.0 

00 

sice 

2.C 

6.0 

0O 

sice 

-1.5 

10.5 

0o 

ccsire 

1.C 

14.0 

450 

This   case  illustrates  what  happens  when  leakage  is  present. 

The  data  window  truncates  the  f/f   =   10.5   component   at   a 

s 

non-integer   multiple   of  the  wavelength,  which  is  the  cause 
of  leakace. 

Case   3->   This   waveform   consisted   of   the   following 
ccn  pcnents 


Amp  . 

f/f 

s 

Phase 
Shift 

E-C 

0.8 

0.0 

00 

cosine 

1.C 

2.0 

0O 

sice 

2.C 

6.0 

00 

cosine 

1.0 

14.0 

450 

ccsire 

-2.C 

22.0 

00 
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This  *a»€fcra  has  a  component  greater  than  (f  )/2,   thus 

s 

this  case  was  used  tc  illustrate  the  aliasing  phenoiencn  and 

the  effect  S  (being  odd  or  even)  had  on  the  folding  cf   this 
component  into  the  analyzing  fcandwidth. 
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APPENDIX  B 
EBIEF  REVIEW  OF  INTEGERS  AND  INDICES 


The  fcllcwing  definition  and  theorems  are  provided  fcr  a 
tetter  understanding  of  the  derivation  of  the  prime 
transfers  alccrithi. 


DEFIIIIICKS 


Onit->  a  unit  is  an  integer  that  divides  every  integer. 
Since  +1  and  -1  divide  every  integer,  they  are  both  units. 

Null  element-*  it  is  an  integer  that  divides  only 
itself.  Kc  integer  different  frcm  zero  is  a  null  element. 
Zero  is  the  null  element  of  the  rational  integers. 

Associates  of  an  integer-*  they  are  the  results  of 
multiplying  the  integer  by  the  units. 

Prime->  it  is  an  integer,  not  a  unit,  that  is  divisible 
by  enly  its  associates  and  the  units.  This  definition 
implies  that  the  greatest  common  diviscr  of  a  prime  "p"  and 
an  integer  "a"  is  1,  or  the  pesitive  associate  of  "p". 

Composite-*  it  is   an   integer   that   is   not  the   null 

element,  a  unit,  cr  a  prime,  e.g.  the  integer  256  is  said  to 

be  highly  cenpesite  while   the   integer   563   is  moderately 
cempesit e. 
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Belatively  prime-*  two  or  more  integers  are  prine  to 
each  ether  cr  relatively  prine,  if  their  greatest  ccamcn 
divisor  is  *1.   The  integers  6,-9,14  are  relatively  prime. 

Tctiect  cf  m->  denoted  0(m)  is  the  number  of  positive 
integers  less  than  "m"  and  relatively  prime  to  "m". 


E.   THEOBEKS 

1.  Iwo  integers  "a"  and  "b"  are  said  to  be  congruent 
modulo  and  integer  "m"  if  "m"  is  an  integer  divisor  cf  a-b. 
This  relaticr  between  "a"  and  "b"  is  incst  often  denoted  by 

a  s  b  mod  (m)  (E-1) 

and   is   read  "a  is  congruent  to  b  mod  (m) "  Eut  for  the  ease 
cf  nctaticn,  it  is  denoted  as 

(  {a)  )  =  b  mod  (m)  (E-2) 

in  the  derivation  cf  the  prime  transform  algorithm. 

If  "a"  is  congruent  to  "b"  modulo  (ra)  and  0<b<mf"b"  is 
called  the  residue  of  a  modulo  (m) .  If  "m"  is  net  an 
integer  divisor  of  a-b,  "a"  and  "b"  are  said  to  be  distinct 
medulo  m. 

2.  If  "a"  and  "a"  (<0)  are  relatively  prime  integers, 
then 


(  {/W  ))  =  1  mod  (m)  (E-3) 
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Corollary  of  2 .   If  »p»   is  a   prime   which   is   not  a 
divisor  of  ai  integer  "a",  then 

F-1 
((a   ))  =  1  mod  (p)  (E-4) 

If  "a"  and  "iM(>0)  are  relatively  prime  integers,  and  Mu"  is 
the  least  positive  integer  for  which 


u 
<(a  ))  =1  mod  (u)  (E-5) 


"uM  is  called  the  exponent  to  which  "a"  belongs  modulo  (m) . 
If  the  exponent  to  which  "a"  belongs  module  (m)  is  #  (m)  ,  "a" 
is  defined  as  a  primitive  roct  modulo  (m)  (or  a  primitive 
icct  cf  b)  . 

3.   The  enly  integers  which  possess  primitve   roots   are 

n    r 
2,4,p  #2p  ,   where   "p"  is  an  odd  prime.   If  "m"  doss  have  a 

primitive  root,  it  has  #(0(m))  distinct  roots  module  (a). 

Let  "pn  te  any  prime,  and  let  "r"  be  any  primitive  root 
cf  "p".  lo  each  integer  "a11  relatively  prime  to  "p"  there 
corresponds  a  unigue  integer  "i"  such  that 


((a))  =  r1  mod  (p)    (0<i<p-1)  (E-6) 


The   integer   "i"   is  called  the  index  to  the  base  "r"  cf  a 
modulo  (p) .   Ihis  is  written  as 


i  =  ind  a  (E-7) 

r 
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Properties  of  Indices 


(a)  ind  1  =  0  (E-8) 


(t)  ind  (-1)  =  <p-1)/2  (E-9) 


(c)  ind  (afc)  s  ind  a  +  ind  b  mod  (p-1)  (E-10) 


D 

(d)  icd  a  s  n  ind  a  mod  <p-1)  (B-11) 


The   power   residue   of  an   integer   "b"  to  the  base  "r"  is 
simply  the  integer  ws"  such  that 


b 
((s))  =  r   mod  (p)  (0<s<p)  (B-12) 


SiDC€  the  power  residue   of   ind  a   is   congruent   to   a 

r 

module       (p) ,      the      concept    of    indicies    and    power  residues   is 
analogous    to    that    of    logarithms    and   antilogrithms. 
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APPENDIX  C 


PRGGRAM  LISTING  OF  THE  FFT  ALGORITHM 


C 

C   THIS  PROGRAM  USES  THE  CEC IMATICN-IN-T IMS  ( CCCL EY-TUK EV  ) 

C   FFT  ALGORITHM  TC  COMPUTE  THE  DISCRETE  FGURIcR  TRANSFORM 

C   (CFT)  COEFFICIENTS  OF  A  COMPLEX  INPIT. 

C 

CIMENSICN  TV(1024) 

C  PENSION  FV(514) 

CIMENSKN  YU024) 

CIMENSICN  AMPQ024)  ,PHASE(6GC) 

CINENSICN  S ( 1012 )«INV (1012)  ,M(2  ) 

CIMENSICN  TIME(60C)t>S(600) 

CCMFLEX  Z 

CCMFLSX*8  A(1024,l,l) 

INTEGER**  ITB(12)/12*G/ 

PEAL*4     RT8(23)/28*0.G/ 

C  REAC  IN  Ti-E  POWER  OF  TWO  "KK"  WHICH  CETERMIN6S  THE  NUM- 
C  ESP  CF  CATA  PCINTS  ANC  CCMPLEX  FRECLENCY  COMPONENTS  ANC 
C   TI-E  LENGTH  OF  THE  CATA  WINCGW  "T"  IN  SECONDS. 

GEAC(5,£  )KK,T 
t       FCFMAT(  I5,E10.7) 
N  1=2**KK 
N2=l 
N3*l 
NMM/2 
IS    *    T/M 


C   CCMPLTE  ThE  SAMPLED  VALUES  OF  AN  ARBITRARY  REAL  WAVEFCPM 
CC  215  J  =  1,M 

XStJ)1"*  0.8  +  C0S(6. 28318  *  2.  *  K  *  TS  )  +  ( 2.0  )*S  IN  (6 
12  E  -  IS  *  6*0  *  K  *  TS  )  + 

S        C^S(<6. 23218  *  14.  *  K  *  TS)  +  Q.7852S) 
$-  (2.0)*CCS(6.28218  *  22.0  *  K  *  TS) 
215  CONTINUE 

M  U)=KK 

* (2  )=0 

M3)  =  0 

C   CAL<"LLATE  THE  FUNDAMENTAL  FREQUENCY  ANC  ITS  HARMONICS. 
C   CALCULATE  THE  TIME  AXIS  VALUES. 


FLNC  '     (l./i  ) 

FV(l)=C.O 

CC  20  J=2tNH 

FN,  ( J)=FV(  j-D  +  FUNC 
2C  CCNTINLE 

TV(1)=0.0 

CC  20  X=2,N1   _ 

TWK)-TV(K-1)+TS 
2C  CCNTINLE 
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C   SET  7FE  CQMFLEX  "A"  MATRIX  TC  2ERC. 


C 

CC    2  11  =  1, M 

CC    2  12*1, N2 

.    CC    2  I3=1,N3 
*    £rnltI2.!3)»(0.0,0.0) 

bJ:     =  iisl?M 

CC    2  12=1, N2 

CC    2  13  =  1, N2 

C       I^NSFER    TFE    SAMPLED    INPUT    DATA    TO    TFE    COMPLEX    "A"    MA- 

cc  3  him**1  =  XS(m 

C   SSfC17CAtJ!  0ISCRETe  FCURIER  TRANSFORM  OF  THE  SAMPL5C 

c  C/LL    FAPM(A,M,INV,S,-1,IFERPI 

C       FclsEMsT5    TH£    VAG,NITUC£    ^C    PHASE    CF    THE    FRSCUENCY    CCM- 

CC    23    11=1, NH 

Z    =    Mil,  1,1) 

ZPEAL    =    REAL(Z) 

ZINAG    =    AIMG(Z) 

I  F(ABSdREAL)  .LT.O.COl)    ZREAL    =    CC 

FFASEdl)    =    1C0.0 

a-   r^SfJ.1!1  =  <ATAN2(ZIMAG, ZREAL))  *  57.2957 
_j   L  C  n  T  I  NL  c 

C   CL7PLT  TFE  COMPUTED  CATA. 

/n    WFI7E(6,47)  T,TS,M,FLNC 

47  FCFMAK/,1  SAMPLE  TIME  =  «,F7.4,'  SECONDS'//, 
i«  S^PLE  INTERVAL  =  «,F6.4  '  sIcCNc!'//* 
5'  NLMEEP  CF  SAMPLES  =  ',14//, 
J'  FLNCAMENTAL  FRECUENCY  =  •,F6«1,»  HZM 
[ T  E (3  )  =  8 
HE  (4)  =  5 
I7E  (12)=1 
NF  =  Nl 
PTE  (1 )  =  C.l 
F7E(2)  =  2.0 

C   FLCT  TFE  SAMPLED  INPUT  WAVEFORM. 

CALL  CPAWP  <NH,7V,XS,I7E,R7E) 

NF  =  Nl/2~ 
PTEdl  =  2.0 
P7E  (2)  =  0.2 

C   FLCT  TFE  MAGNITUDE  OF  THE  CNE-SIDEC  SPECTRA. 


CALL  CRAWP(NH,FV,AMF  ,IT6,RTE) 
P7E(2)  =  90.0 


C   FLCT  TFE  FFASE  GF  THE  CNE-SIDEC  SPECTRA. 

CALL  CRAHP(NH, FV, FFASE, ITE, RTB  ) 

S7CF 

ENC 
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APPENCIX  C 


PROGRAM  LISTING  OF  THE  CZT  ALGORITHM- 


?cIfi,F5SS?5K»£2WEyiES  THE  C*SCPETE  FREQUENCY  CCMFCNENTS 

^   c    I7a    L52fLE)(  INPUT  *AVEFCW  BY  THE  CHIR°-Z- 
IrANSFCR^  ALGCPITHM. 

XJig24),Y<1024),XW(512),YW<512>,XS(512) 

YS  512,XP<512),YP<512) ,FREC(512),WR<512) 
£Mfi?h^(fi?MI(512),AMF(512),TIME<512) 
fHASE<512),SS(512),F,T,N,FN,SPL,PRPW,TRW,PCRW 
FCT1,FCT2,FCT3,FCT4,FCT5,FCT6 

ITB(12)/12*C/ 

PTe<28)/28*C.C/ 


CCNNCN 

CC*MGN 

CCNNCN 

CCN^CN 

CCMNCN 

INTEGEP*4 

PEAL*  4 


FEU  THE  PARAMETER  "II",  WHICH  CETEPKINES  IF  CATA  IS  TO 
EE  FE*C  IN  CR  CALCULATED  WITHIN  THE  FPCGRAM.   IF  "II"  IS 
ECLAL  TC  "C",  THE  PRCGRAM  WILL  GENERATE  A  CCNPLEX  SINL- 
SCICAL  INFLT  FCR  A  SINGLE  FREQUENCY.   IF  "II"  IS  NGT 
ECLAL  TC  ZERO,  AN  ARBITRARY  REAL  INFLT  MAY  EE  FEAC  IN  ANC 
T  h  E  ^^G  IN  A  FY  PART  (YS)  IS  SET  TO  ZEPG. 


IK 


IFLN    =    1 
FEAC<5, 110) 
FCFMTH5) 


II 


"IV"  IS  THE  NUMBER  OF  TIME  AND  FREQLENCY  SAMPLE  PCINTS. 
"F"  IS  THE  FRECUENCY  TC  eE  USEC  IN  SLEROUTINE  WAVE  ANC 
■T"  IS  THE  LENGTH  CF  THE  INPUT  CATA  WINDOW  IN  SECCNCS. 
"SFL11  CETEFVINES  THE  PATE  CF  SFIRALING  CF  THE  CCNTCUR 
IN  THE  Z-PLANE.   IF  "SPL"  IS  ECUAL  TC  ZERO,  THE  CCNTCLR 
IS  A  LMT  CIRCLE  IN  THE  Z-FLANE  WHICH  CCRRESFCNCS  TC 
CALCLLATING  THE  DISCRETE  FCURISR  TRANSFORM, 
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FEAC(5,  105)  F,T,SPL,N 
FCFNAT<2E10.7,I5) 


CEAC  IN  THE  PARAMETER  "FRRW"  WHICH  CETERMINES  WHETHER  CR 
NCT  FANDC^  VALUED  FR  E-NLLT IPLY  INC-  CCEFFICIENTS  WILL  EE 
LSEC  IN  THE  CCMPUTATICN.   PARAMETER  "TRW"  DETEFMINES 
AETHER  CR  NCT  RANDOM  VALUED  TAPPING  WEIGHTS  IN  THE  CCN- 
VCLLTICN  FILTERS  WILL  EE  LSEC  IN  THE  CCMPUTATICN. 
F4RANETER  "PORW"  DETERMNES  WHETHER  CR  NCT  RANCCN  VALLEC 
FCST-NLLTIFLYING  WEIGHTS  WILL  eE  USEC  IN  THE  CCVFUTATICN. 

F*AU5,101)  PRRW,TRW,PCPW 
1C1  FCFNAT(2E10.7) 

"FCT1"  AND  "FCT2"  DETERMINE  THE  RANGE  OF  VALLES  THE 

FPE-NLLTIFLYING  COEFFICIENTS  CAN  CEVIATE  FRCM  I    NOMINAL 

WLLE.   "PCT3"  AND  "FCT4"  DETERMINE  THE  RANGE  CF  VALLES 

Tl-c  FILTER  TAPPING  WEIGHTS  CAN  CEVIATE  FROM  A  NCMNAL 

WLLF.   "PCT5"  AND  "FCT6"  DETERMINE  THE  RANGE  CF  VALLES 

FCWLLTIPLYING  COEFFICIENTS  CAN  DEVIATE  FRCM  A  NOMINAL 
V  ALL  E. 
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1C2    PF^J?tJi?JJSTl'     FCT2>    PCT2>    PCT^    PCT5'    FCT£ 
FFEC(l)    =    CO 

Ff>   =  n 

ill  Wli  I  iB  few'BI^'" E£UAL  7C  ZER0 

cc  «o   M  =  1,N 

*'  ■  P»  4  1 

^lf)  =  C.C 
FFEC(J)  *  N  *  (1./7) 
4C  CCMINLE 
GC  TG  S 

r       cjSIcS-II  It!  V$tWIS  0F  *N  ARBITRARY  CCMPLEX  WAVEFCRM 
L   i^VFLcL  "N"  TINES. 

5    C/LL  WAVE 
9    C4LL  FRMLT 
C / L  L  XSLM 
C/LL  TAFW 
CUL  CCNV 
C/LL  FCNULT 
C/LL  ANFhA 
C 
C   CLTPLT  ThE  CCNPUTED  CATA. 

WPITE<6,220)  IRUN,FCT2 
220  FCFNAT</ ,2X, 'THIS  PUN  IS  NUMBER   SI2, 
$♦   ANC  ThE  PERCENT  CEVIATIGN  IS  «,F2.2) 
I7E (2  )  =  € 

11EIM     =    5 

ne  (12)  =  i 
RTE(i)    =  o.i 
PTE (2 )    *    2.0 

KfFTS    -    N 

C 

C   FLCT  THE  SAMPLED  INPUT  WAVEFORM. 

C 

C/LL  CP4WF(NUMPTS,TIMEf >Sf ITEtRTE) 

R1E(2)  =  0.4 

C   FLCT  ThE  IMPULSE  RESPCNSE  CF  ThE  CCSINE  FILTER. 
C 

C/LL  CR/WF(NLMPTS»TIVE,X,ITB,RTE) 
C 

C   FLCT  ThE  IMPULSE  RESPCNSE  CF  ThE  SINE  FILTER. 
C 

C/IL  CFAWF(NUMPTS,TiyE,Y,ITE,R7E) 

iUMFTS  =  N/2 

F7E<1>  =  2.0 

FIE  (2)  *  0.2 

C   FLCT  ThE  MAGNITUDE  CF  THE  FREQUENCY  CCNPCNENTS. 

C/LL    CR/WP(NUMPTS,FPEC,AMP,ITB,PTE) 
FT  E  (2  I    *    SCO 

C       FLCT    7FE    PhASE    CF    ThE    FREQUENCY    CONFCNENTS. 


C 


C/LL    CRAWP(NUMPTS,FPEC,FHASE,ITe,PTB) 

STCF 

EK 
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SIEPCLTINE    WAVE 

c 

C       1HS    SLERCLTINE    COMPUTES    THE    VALUES    OF    AN    ARBITRARY 

C       SELF-GENERATEC    COMPLEX    WAVEFORM    FOR    A    SINGLE    FFECUENCY. 

C       THE     FFSCUEfvCY    IS    SPECIFIEC    BY    THE    LSEP. 

C 

CCfNCN    >{1024),Y(1C24),XW(512),YM512),XS(512) 
COMGN    YS(512),XP(512),YP(512),FREQ(512)  ,WR(512) 
CCNNCN    V»I(512J,GR(512),GI(512),AMF(512),TIME(512) 
CCfMCN    FHASE(512) ,$$(512) ,F,T,N,RN,SPL,FRRW ,TPW,FCRW 
CCf-NCN    FCTl,PCT2,PCT3,PCT4,PCT5tPCT6 

15  *  m 

CC    2C    *  =    1,N 

K    -    M    -  1 

J    =    M    ♦  1 

FF  EC< J)  =    M    *    (l./T) 

XS(M    =  CCS(6. 28318    *    F    *    K    *    TS  ) 

VMM    =  SIM6. 28318    *    F    *    K    *    TS  ) 
30    CCMIME 
FHLFN 
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SIEPCUTINE  PRMULT 

|   8fM™T,M  ^WB  IB  KftMMNAffi/Mi&tt! 

CCNNCN  X(1024),Y(lQ24)tXM512).YW(cl?l.x<;f  =  i:M 

K  »  M  -  1 
C   CCMF11E  7hE  TIME  INCREMENTS  FOP  PLCTTING  THE  CLTFUT. 

TINE(M)  =  K  *  (T/RN) 
C   CCHFITE  THE  PRE-MULTI FLICATIQN  CCE F FIC  IENTS . 

S  »  (SPU*  (K**2)/2. 
CTF  =  E>P(S) 

c  VMM    =    (CTR)    *    5IN<<3.1415?MK**2M/RNI 

?       JJ.'5.*.lliSj    Ifcl       PR5-MLLTIFLICATI0N    COEFFICIENTS    ARE 
£       ?ffcC^JI|£TfQg    ScUNJfCGNLY    RANCCM    <?ANGE    CETEFMNEC    eY 

C       ThE    FcFcfcNTAGS    3  F    DEVI4TICN    FRCP    A    NCMINAL    VALIE. 

IF(FF.RW.NE.l.O)    GC    TO    10 

CXLL    P/UCUIX,IY,YFl) 

PhC    =    VFL 

>MM    =    Xk(H)    -<  (PCT1)*XW<M )  )    +     (PNC*PCT2*XMm 

VV(M)     =    YW(N)    -<<PCTI)*YW(M))    +     ( RNC*DCT2* Y M V  )  ) 

1C       CCMIME 
FETIFN 
EK 
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SLEPOLTINE  XSUM 
C 

C  ThIS  SLEPCLTINE  MULTIPLIES  ThE  SAMPLEC  CATA  ev  ThE  PRE- 
C  NLLTIFLICA7ICN  COEFFICIENTS  ANC  SU^S  THE  R5*L  4NC  I^AG- 
C   IN/RV  COMFCNENTS. 

CCNNON 

ccn*cn 

CCNNCN 

CCNKN 

CC  16  M  =  1,N 

>F1  =  XMM)  *  XS(M 

>F2  =  XS(M)  *  (-1.  *YWtMU 

Yll  =  YS(M  *  (-1.  *YMM)> 

YI2  =  YS(M)  *Xfc(V) 

C   3LM  ThE  REAL  AND  IMAGINARY  TERNS  CF  ThE  PRE-NL LT I  PL  I- 
C   C4TICN  PPICP  TC  CONVOLVING  IN  ThE  ChIPP  FILTEPS. 


C 


fcMM    '    XR1    -    YU 
H  (N)    =    XR2    +    YI2 
16       CCMINLE 
FE7LRN 
EK 
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S  L6PCLT  INE  TAPW 

TUS  SLERCLTINE  COMPUTES  THE  TAPPING  WEIGHTS  FCP  THE  CCK- 
FLEX  CCNVCLLTICN  FILTERS. 


11 


14 


CCNNCN 
CCNNCN 
CCNKN 
CCNKN 
CCPM2N 

pr   *  n 

CC    11   f    =    1,N 

K    -    N-    N 

$    =  <(SFU*    (K**2J/2.) 

C1P    =    EXP(S) 

X(M       =     (CTR) 

V  <>  )       =    (CTR) 

ccminle 

NM    =  N  -    1 

CC    14  K  =    1.NM1 

J    »   N  -  K 

>  -    N  ♦  K 

>  <M    =    >(J) 

>  < N  )    -    Y(J) 
CCMINLE 


1(512),  YM512)fXS(512) 
»(512),FBEC<512)  ,UP(512) 
<512),AMP(512) ,TIME(512) 


X(1024),Y(1024),XW 

YS(512),XP(512),YP, 

WI  (512),  GR  (512),  GHdiij  v  Pi"  P13J.ZI  »  i  imti  sn  j 

FHASE(512),SS(512),F,T,N,RN,SPL,PRPV»,TPW,FCP* 

FCT1,PCT2,FCT3,PCT4,FCT5,FCT6 


*  (-1.0) 


*  CCS((2.1415S*(K**2))/RN) 

*  SIN((3.14159*(K**2)  )/RN  ) 


IF  CAllEC,  TFE  TAPPING  WEIGHTS  ARE  CALCULATEC  FCP  A  UM; 
FCFK>  RANCCN  OANGE  CETSRNINEC  BY  THE  PERCENTAGE  CF  CEVI- 
A  T IC  h  FPCN  A  NCMINAL  VALUE. 


IF( 

N  P  W 
IX 

cc 

C  A  L 

pre 
>(>J 

>  (J 
I> 

CO 
FET 
EK 


TP  W 

■*1 

15 

I  P 
1  * 
)  * 
)  = 

■  I 
TIN 
LRN 


.NE.l 

2  *  N 

11111 

J  =  1 

ANCL( 

VFL 

>(J) 

XJ) 

LE 


.OJRETLPN 

-  1 
111 
,NRW 
IX,IY,YFL) 

-  ( (PCT2)*X(J) ) 

-  ((PCT2)*Y(J)) 


+  (BNC1*PCT4*>( J)  ) 

+  (pm:i*pct4*y(  j)  ) 


168 


SlEP-CUTINE  CCNV 
J£2?rTi"STi=tTINE  SHIFTS,  MULTIPLYS,  ANC  SUMS  7C  FERFQRN 

ccn vc ilticn. 

uf  IS,1?  '££  IH  \'XPr  <512),FREC(5i2l.*»R<512) 


cc 

CC 
CC 
CC 
CC 
CC 
SL 
SI 
SI 
SI 
CC 
IF 

Gl 

G2 
G2 
G< 
GC 

Gl 
C2 
G3 
G4 


MCN 

nngn 

15  N 

M  = 
M  - 
M  - 
M  - 

in   k 

(N.GT 

»  N  ♦ 


TC  3 

N  + 


CO 
CO 
CO 

CO 
=  1,N 

•  1)  GG 
1  -  K 
WP(K) 
U=  (K) 
U!(K) 
WI(K) 

v  -  K 
W  P  (  K  ) 
UR<K) 
WI(K) 
WI(K) 


TO 


X(J) 
Y(J) 

X(J) 
Y(J) 


X(J) 
Y(J) 
X(J) 
Y(J) 


SUN  LF  THE  VALUES  AFTER  SHIFTING  ANC  MULTIPLYING  IN  THE 
CCNVCUTICN  FILTERS. 


SIM  =  Gl 
SLN2  =  G2 
SIM  =  Q2 
SL^^i  =  G4 
CCMINUE 


SUM 

SUM 
SUM 
SUN  A 


CCVFL1E  THE  REAL  AND  IMAGINARY  CCMCNENT  VALUES  AT  THE 
ClIFU  CF  THE  CGNVCLVEP. 

GF  (M  =  SUM  -  SUM 
Qlif  )  =  SLM  +  SUM 
If   CCMINUE 
FE1LPN 
EK 
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c      51EFCLT1NE  PCMULT 

r   Jn^fHf^  Itf  POSJ-fLLTIFLICATIGN  CCEFFICIEM*  ARE 

CfU    R^CU(IX,IY,YFL) 
FMI*    s    >  FL 

I>    S    IY  -<(PCT5)*YP(MI)    ♦     (PND2*PCT6=»VF(?)) 

65       CCf^TIME 
PETtRN 

Ere 
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SIEFCLTINE  AMPHA 
C 

C   THIS  SLBPCLTINE  PCST-NLLTI  FLI ES  THE  REAL  ANC  UAGINAPY 
C   CCNVCILTICN  FILTER  OUTPUTS  AND  THEN  CETERMINES  THE  PAGM- 
C   TLCE  ANC  PHASE  CF  THE  FPECUENCY  CCI*FCNENTS. 

CC*KN  X(1024),Y(1024),XH(512),YM512),XS<512) 

CO  VON  YS(512),XP(512),YP(512),FFEC(512),V»P<512) 

CCNVCN  WI(512),GR(512),GI(512),ANF(512),TIVE(512) 

CCNKN  FH  AS  E  <  51 2  ),SS(  5  12  ),F,T,N,PN,SPL,PPPV«, TRW,  FCRW 

CO  VON  FCTl»PCT2tPCT3,PCT4, FCT5t FCT6 

CC    SO    f    =    1,N 

ei    *    GF(N)    *XW(M> 

£2    *    GR(P)    *    (-1.    *YW(M)) 

£2    =    GUM    *    (-1,    *YW(M)) 

E«    =    GI  (I*  )    *XW(P) 

C    SLf    TFE    FCS7VLLTIPLIEC    VALUES    TC    CETAIN    THE     IMAGINARY    ANC 
C       PEAL    VALUES    CF    THE    FPECUENCY    COPCNENTS. 


C 


PE  =    ei    -    83 

PI  =    E2    ♦    B4 


C   CETEFMNE  THE  MAGNITUDE  ANC  PHASE  CF  THE  FREQUENCY  CCPPC- 
C   NBMS. 

AFC-    =    (PE**2)    +    ( P  1**2  > 

^F(M)     *     {(SGRT(ARG)  )/RN) 
IF(AES(PE).LT.O.OG1)RE    =    O.C 
IF  (AES(FI).LT.O.OOliRI    =    0.0 
FFASE(M)    *    100.0  „    „ 

IF(FE.EC.0.0.ANC.RI.GE.C.0)PHASE(N)    =    90.0 
IF (RE.EC.G.Q.AND.PI.LT.O.O)PHASE(M    =    -SCO 
IF(PE.GE.G.O.ANO.RI.EQ.C.Q)PHASE<M    =    0.0 
IF  (FE.LT.O.O.AND.FI.EC.C.O)FHASE(f«  )    =    18C.0 
IF(FHASE(N).NE. 100.0)    GC    TC    SO 
FFASE(M)    =    ATAN2(PI,RE)*    57.2S57 
SC       CCM1NLE 
PE1LPN 
EN  C 
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APPENCIX  E 


FRCGPAN  LISTING  OF  THE  PRIME  TRANSFGRM  ALGCRI7FP 


ThIS  FROGRJM  USES  THE  "PRIME  7PANSFCPM"  ALGORITHM  TC 
CCMFL76  TFE  FCURIER  COEFFICIENTS  OF  SAMPLED  PEAL  CATA. 

CCNNCN  >S(32),X(64),Y(64),P(32,32),FT<32,32),GP<32) 
CCNKN  NREM(32},Rc<32),PI(32)tXMAC-(32),X?<32) 
CCNNGN  FRE<32),PRI<32),FHASE(32),TIME<32),FFEC<32) 

ccm^cn  n,t,nn,ts 

IhTEGEP*4  ITE(12)/12*0/ 
FEAL*4     RTB(23)/28*O.Q/ 

FEAC  IN  THE  CATA  WINCCW  LENGTH  "T"  ANC  THE  NLVEER  OF  CATA 
SAfFLES  "N",  WHERE  "N"  IS  AN  GCC  FRINE. 

FEAC(5,10)T,N 

FEAC  IN  THE  RtSICUE  V&LLE3  TAKEN  <=RCN  CRC  TABLES  FOR 
HE  FFLATICN  (R**N)  KC(N'). 

^^  =  n  -  l 

PE AC ( 5,  11)  (NREM(K),K  =  1,NM) 

CALL  S'aAVE 

CALL  FRNLT 

CALL  PEFfcAV 

CALL  TRANU 

CALL  TPCSE 

CALL  )i^QP\- 

CLTFLT  THE  CCNPUTEO  CATA. 

WFITE (6 t205) 
WFITE(6  ,24C) 

WFITE(cUoO)((P(  I»J)»J=ltNM)  ,I=1,NM) 

U  FI T6     (6,240) 

fcFI7E(6,215) 

£fItI(6 jloo!(<PT( I,J) »J«1,NH1  ,I«1,NM) 

fcFI7E(6,24C) 
WFITE(6,260) 
*FI/7E<6,22G) 

SfIT6(6;Iio!(I,XS(I),XP(I),I=1,N) 

WFI7EU,2A0) 
WFI7E(6,235) 
WFITE (6  ,245) 

tpiTi(6!230)(I,NREM(I),X(I),>r{n,PE(I),RI(IJ,FRE(I), 

5FCI  (I  ),I  =  1<,NM) 
*FI7516,24C) 
WFI7E(c,255) 

kFITEUfllOM  I,XMAG(I),PHASE(I)  ,I=1,N) 
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ITEt2)  =  2 
ITE(2)  *  8 
I7E(4)  *  5 
I7E  (12)  =  C 
PTfU)  =  T/8.0 
R7E(2)  =  2.0 

FLCT  THE  SAMPLED  INPIT  WAVEFORM. 

C^LL  CPAWP(N,TIM6,XS,ITe,RTe) 
F7E  (1 )  *  2.0 
NFT  =  N/2 
RTE(2)  =  0.2 

FLCT  THE  MAGNITUDE  OF  ThE  FREQUENCY  COMPONENTS. 

CUL  CR4WPCNPT,FRECtXMAG,ITE,RTB) 
**  T  c  1 2  J  s  Su.O 

FLCT  TFE  PHASE  CF  THE  FRECLENCY  CCNFCNENTS. 


1C 

11 

110 

2CC 

2C5 

21C 

215 

22C 

23C 


<C 


25  £ 

2tC 


CAL 

FCF 
FCF 
FCF 
FCF 
FCF 
FCF 
FCF 
FCF 
FCF 
FCF 
$'SI 
$'F5 
FCF 
FCF 
$4>, 
FC  = 
FCF 
STC 
ENC 


b  c?^i;(NpT» prophase, iTe,RTE) 

NAT( 510.5 tI5  ) 

*AT{  1615) 

MT(/,I5,2F20.2) 

NAT(/,2X,30F4.1) 

MATC40X, 'PERMUTATION  MATRIX  «F«") 

NAT(/,2X,I5,3X,2F20.5) 

MAT(25>,«THE    TRANSPOSE    OF    THE     FERMUTATICN    NATRIXM 

M?i5?S;i?S!rii}i?i^§?ML    """0..13X..PERJ.0TEC.I 

MAT(2X, 'INDEX*  ,    5X , • RESI CUE*  ,5X ,  «CCS INE    WEIGHT • ,5X, 

NE    rtclL-HT*  ,7X,  'REAL    PARTS    6X ,  ♦IMAGINARY    PARTff4Xt 

£L    PART', 6Xf •IMAGINARY    PART") 

*MT(//) 

MAT(57X,«PERMUTEC  CRCER •  ,4X  ,' PERMUTED  CRCEF', 

•NORMAL  CRDER1 ,4X,«NCRMAL  CRCEF1) 

MATC 2X, • INDEX • ,10X,« MAGNITUDE  •  ,12X,' PHASE' ) 

VATC24X,'  INPUT  CATA  •  ) 

F 
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SIEFCLTINE    SWAVS 


TUS    S 
TRAR\ 
THE    SA 
KMCS 

7FE    S^5 
IN    ThE 


CCN 
CCN 
CCf< 
CCN 
7$ 

cc 

K  = 
TIN 
FFE 
XS( 

1262 

$ 

*-  ( 
CCN 
PET 
ENC 


L3SCLTINE  COMPUTES  ThE  SAMPLE  VALUES  FCR  AN  ARBI- 
fcAVEFCRM  AS  DEFINEC  EY  ThE  USER.   IT  ALSC  COMPUTES 
MPLE  TI*E,  THE  FLNCAMENTAL  FREQUENCY  ANC  ITS  FAR- 


NPLE  TINES  AND 
CLTFUT  PLOTS. 


FRcCUENCIES    A*E    USEC    AS    The    X-AXIS 


12 


NCN 
NCN 
NCN 
N3N 
*    T/ 
12    N 
M    - 
E(M) 
C(M) 
M    = 
18    * 

2.0) 
TINL 

LRN 


XS( 

NRE 
PRE 
NtT 
N 


32), X (64 ), Y ( 64 ), P (32, 2 2  ), FT( 3 2, 3  2  ),GP( 22) 
V(32),RE(32),RI(32),XNAG(32),XF(32) 
(32),PRI(32),PHASE(32),TIME(32),F«EC(32) 
♦NN,TS 


=    1,N 


0. 
6. 

CCS 
*CC 
E 


K    *    TS 

K    *    (l./T) 

8    +    C0S(6. 28318    *    2.    *    K    *    TS  )    + ( 2.C ) *S  IN ( 6 

0    *    K    *    TS)    ♦ 

((6.28318    *    14.    *    K    *    TS)    +    0.78529) 

S(6. 28318    *    22.0    *    K    *    TS) 
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SLEFCOTINE  PRMUT 

C   THIS  SieRCLTINE  CGMPUTES  THE  PERMUTATION  MATPI*  "P"  ThftT 
C   IS  LS6C  TC  PERMUTE  Tl-E  SAMPLED  REAL  CATA.   IKE  TPANSFCSE 
C   Cf  "P"  IS  ALSC  CGMPLTEC. 
C 

CCNNGN  XS(32),X(64),Y(64),  d  (  32  ,  22  )  ,PT  (  32  ,  22  » ,  C-P  (  2  2  ) 

CCNKN  NREM(32),RE(22I,RI(32),XNAC-(32>»XF(22> 

C CM  NGN  FRE(22),PRI(22),PHASE(22>,TIME(32),FPEC(32) 

CCNfON  N,T,NV,TS 

CC  20  NN  *  1,NM 

CC  25  M.  =  1,NM 

IF(NPEM(NN).N5.M)  GC  TG  1 

F(NN,M)  =  1.0 

GC  TC  25 
1    F(NfwM)  =  CO 
25   CCMINLE 
2C   CCNTINLE 

C   TPANSPCSE  THE  PERMUTEC  MATRIX 
C 

CC  E5  M  =  1,NM 
CC  <0  K  =  1  ,NM 
FT(K,N)  =  P(M,K) 
<C  CCMINLE 
Ef  CCMINLE 
PETLPN 
ENC 
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SIEFCL7INE  PERWAV 


HIS  SlERClTlNE 

FIRST  SANFLE. 

JtLFilS  IC  6E  PERMUTEC  IS  AN 

E!iriTf?  ?V  MULTIPLYING  17  BY 

^iCt  H  AN  (NM  X  W>    matrix. 
Iff    X  1 )  • 


PERM7ES  THE  SAVPLEC  CATA  EXCUCING  THE 


(NM  X  1)  MATRIX  *NC  IS 
7h=  PERM'TATICN  ^ATRIX 
THE  RESULTAN7  MTIX  IS 


NR^i?:^('||!^}i^^,32,iPT(32,^,,C.P(32, 

fi^zj.fRKaaiJPHAiiJlIj^^iJliJJNlcoz, 

1,NM 
1,NM 


CCNNCN 
CCNVCN 
CCNNON 

ccpkn 

CC  50  I 

Gl  *  C.C 

CC  £0  J  = 

K  =  J  +  1 

PC  =  P(I,J)  *  XS(K) 

G 1  =  P  G  +  G 1 
tC    CCN7INUE 

GP  (II  =  Gl 
50  CCN7INLE 

CC  30  J  =1,NM 

K  =  J  +  1 

>F(1)  =  XS(1) 
„„  XF<K)  ^  GP(J) 
rC       CCMINLE 

FE71PN 

Ef^C 
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SLERCUTINE  TRANW 

TUS  SLSROLTINE  COMPUTES  THE  TAFPING  HEIGHTS  FCP  THE 
CCMPLE*  CIFCLLAR  CORRELATION  THAT  y/!Y  BE  IMPLEMENTED 
L  SING  fi    TRANSVERSAL  FILTER  AND  CCMFLTES  THE  RESLLTS 
CP  TFE  CORRELATION.   THE  TAPPING  WEIGHTS  CCRRESPCNC  TC 
THE  VALUES  OF  THE  »C"  MATRIX. 


CCMMON 
CCMMGN 
CCMNON 
CCVMCN 


>S(32),X(64),Y(64),P(32,32)tPT(32,22),GP<32) 
NREV(32),RE(32),RI(32I»XVAG(32),XF(32) 
FRE(32),PRI<32),FHASE(32),7IME(32),FFEC<32) 
N,T,NM,TS 


C4LCLLATE  THE  TRANSVERSAL  FILTER  COEFFICIENTS. 

CC  10  K  =  1,NM 
I  -    NM  4  K 

X(K)  =  CCS((6. 28313  *  NPEM(K))/N) 
Y  (K  »  =  (-1)  *  SINU6. 28219  *  NREHKM/N) 
XI)  -  >(K) 
YU)  =  Y(K) 
7C   CCNTINLE 

FERMLTED  DATA  IS  CORRELATED  IN  THE  TRANSVERSAL  FILTER. 

CC    7  5    M    =    1,NM 

SLM    =    CO 

SLM2    ~    CC 

CC    60    KK    =    1,NM 

J=NM4l+M-KK 

XF 1  =  X(J)  *  GPCNM+1-KK) 

VF1  =  Y(J)  *  GP(NM+1-KK) 

SLN  AFTER  SHIFTING  ANC  MULTIPLYING 

SLM1  =  >R1  4  SUM1 
SLM2  =  YR1  4  SUM2 
EC  CCMINLE 

THE  FEAL  PA«T  IS  OBTAINED  FROM  THE  CCSINE  FILTER  ANC  THE 
IMGIMFY  P4RT  IS  OBTAINEC  FROM  THE  SINE  FILTEP. 

PE  (M  -  SLM 
P  KM)  =  SLV2 
15  CCNTINLE 
PETLRN 
END 
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SIEFOLTINE  T-POSE 
C 

C   1HS  SLBFOLTINE  MULTIPLIES  THE  CIRCLLAR  CORRELATOR  RE- 
C   SILTS  EY  THE  TRANSPOSE  CF  THE  PERMLTATION  MATRIX  ANC 
C   THEN  ACDS  THE  FIRST  CATA  SAMPLE  TG  EACH  COMFlTcC  REAL 
C   F*P7. 

CCNM3N  XSC32l1X(64lffY(64lfP(32ff32ltPTC32f22It6P(32J 
CChKN  NREM<22),RE(22),RI(32>,XMAC(22>,XF<22) 
CCN^CN  FRE(32),PRK32),FHASE(22),TIME(32JtFREC(22) 
CCNNON  ^T,W,TS 
CC  «5  t*    =  ItNM 

C-l  =  c,c 
GZ  =  O.C 

CC  100  K  =  1,NM 

XI  =  PT(V,K)  *  RECK) 

>1  *  PT<y,K)  *  RI(K) 

CI  *  XI  ♦  Gl 

Q2    =  Yl  ♦  G2 
ICC  CCMINLE 

FFE(M)  *  Gl  +  XS(1) 

FFKMI     *    G2 
SI       CCMINLE 

PETLBN 

Ef^C 
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SIEFCITINE  XMGPH 

C  THIS  SUBROUTINE  COMPUTES  THE  MAGNITICE  AND  PHASE  FOR 
C  5*0  FREQUENCY  COMPONENT  TC  GIVE  A  CI^E-SIDEC  SPECTRA 
C   NCTE  THAT  THE  C-C  COMPONENT  IS  COMFLTEC  SEPARATELY. 


C 


CCfNGN    XS(22),X(64),Y(64),P(32,22),PT<32,22),GP<22) 

CC^MON    NREM<22)tRE(22)tPI(32)tXMAG(22)»XF(22) 

CCfNCN    PRE<22),PRI(22),PHASE(22),7IME(32)»FREC<32) 

CCNNCN    h,T,NN,TS 

CC    115    *    =    1,NM 

J  =  NU 

IF(*ES<FRE(M)).LT. C.OOl)    PRE(M)    =0.0 

IF(4BS(PRI<M)). IT. 0.001)    PRKMJ    =    CO 

*FG    =    (FRE(M)**2)    +    <PRI(M)**2) 

>^4G( J)    =    <SCRT(ARG))/N 

FMSE(J)    =    1C0.0  „    m 

IF(RRE(M).EG.O.Q.ANC.PRI(M).GE.C0)PHASE(J)    =    SCO 

IF(FRE(M).EQ.O.O.AND.PRI(M).LT.CC)PHASE(J)    =-90.0 

1F(FRE(M).GE.O.O.AND.PRI(M).EO.O.O)PHASE(J)    =    CO    m 

IF  (FRE(M).LT.0.Q.ANC.PRI<M).EC.CC)PHASE(J)    =    180.0 

IF<FFASE(J).NE. 100.0)    GC    TO    115 

)    =    ATAN2(PRI(M),PRE(N)  )    *    57.2957 

AND    PHASE. 


c 
c 
c 

115    CCNTINUE 

CCMFITE    THE    "C-C"    COMPONENT    MAGNITUCE 

FFASE(l)    =    CO 

GC    =    CC 

CC    105    w    =    1,N 

GC    s    XS  (Jl    +    GO 

1C5    CCI^TIME 

>MAG(  1)    =    GC/N 

FEUPN 

EK 
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